Problem with ordinary differential equation

Hello, can anyone clarify me why I'm getting a line chart with the following code and how can I fix that? My initial conditions are y(C0)=t0, and y'(C0)=0. Thanks

syms y(x) x Y

N=5;

r=0.05;

m=0.01;

p=1;

s=1;

t=0.1;

a=0.25;

b=0.25;

C0=1;

C=5;

t0= (((1-a)/a)*N^(a+b))/(r-m);

Dy = diff(y);

D2y = diff(y,2);

ode = y-((((1-a)/a)*N^(a+b))-x*m*Dy+0.5*D2y*(x^2)*(s^2))/(r-m);

[VF,Subs] = odeToVectorField(ode);

odefcn = matlabFunction(VF, 'Vars',{x,Y});

tspan = [C0 80];

ic = [t0 0];

[x1,y1] = ode45(odefcn, tspan, ic);

figure

plot(x1,y1(:,1), 'DisplayName','(x_1,y_1_1)')

grid

Answer 1

[Star Strider]

The initial condition of 0 for keeps it 0 througuout the integration.

syms y(x) x Y 
N=5;
r=0.05;
m=0.01;
p=1;
s=1;
t=0.1;
a=0.25;
b=0.25;
C0=1;
C=5;
t0= (((1-a)/a)*N^(a+b))/(r-m)

t0 = 167.7051

Dy = diff(y);
D2y = diff(y,2);
ode = y-((((1-a)/a)*N^(a+b))-x*m*Dy+0.5*D2y*(x^2)*(s^2))/(r-m);
[VF,Subs] = odeToVectorField(ode)

VF = 

Subs = 

odefcn = matlabFunction(VF, 'Vars',{x,Y})

odefcn = function_handle with value:
    @(x,Y)[Y(2);(1.0./x.^2.*(sqrt(5.0).*-3.0e+2+x.*Y(2)+Y(1).*4.0))./5.0e+1]

tspan = [C0 80];                                     
ic = [t0 0];                                     
[x1,y1] = ode45(odefcn, tspan, ic);
 
figure
plot(x1,y1(:,1), 'DisplayName','(x_1,y_1_1)')
grid

.

[Fatemeh]

Thanks much. Instead of defining an initial condition for y2, can I define two conditions for y1? like y(C0)=t0, and y(80)=0. How can I define two conditions for y1? Thanks again

[Torsten]

Your problem becomes a boundary value problem instead of an initial value problem.

Instead of ode45, you will have to use bvp4c.

syms y(x) x Y 
N=5;
r=0.05;
m=0.01;
p=1;
s=1;
t=0.1;
a=0.25;
b=0.25;
C0=1;
C=5;
t0= (((1-a)/a)*N^(a+b))/(r-m);
Dy = diff(y);
D2y = diff(y,2);
ode = y-((((1-a)/a)*N^(a+b))-x*m*Dy+0.5*D2y*(x^2)*(s^2))/(r-m);
[VF,Subs] = odeToVectorField(ode);
odefcn = matlabFunction(VF, 'Vars',{x,Y});
bcfcn = @(ya,yb)[ya(1)-t0;yb(1)];
xmesh = linspace(C0,80,150);
solinit = bvpinit(xmesh, [0 0]);
sol = bvp4c(odefcn,bcfcn,solinit);
plot(sol.x,sol.y(1,:))

[Star Strider]

@Fatemeh — As always, my pleasure!

@Torsten — Thank you!

(I just now saw this.)

[Fatemeh]

Thanks so much!