Need to solve this with Runge-Kutta 4th order method in MatLab.

y''-μ(1-y^2)*y'+y=0

y(x=0)=1 ; y'(x=0)=0

μ=0; 0.2; 1; 5 and 10 (for each value a different result)

y'(x) versus y(x) and x versus y(x) functions plot in the range [0; 50].

[James Tursa]

What have you done so far? What specific problems are you having with your code?

[Cubii4]

fy=@(x,y,z) z;

fz=@(x,y,z) (1-y^2)*z-y;

x(1)=0;

z(1)=0;

y(1)=1;

h=1;

xfinal=10;

N=ceil((xfinal-x(1))/h);

for j=1:N:xfinal

x(j+h)=x(j)+h;

k1y=fy(x(j),y(j),z(j));

k1z=fz(x(j),y(j),z(j));

k2y=fy(x(j)+h/2,y(j)+h/2*k1y,z(j)+h/2*k1z);

k2z=fz(x(j)+h/2,y(j)+h/2*k1y,z(j)+h/2*k1z);

k3y=fy(x(j)+h/2,y(j)+h/2*k2y,z(j)+h/2*k2z);

k3z=fz(x(j)+h/2,y(j)+h/2*k2y,z(j)+h/2*k2z);

k4y=fy(x(j)+h,y(j)+h*k3y,z(j)+h*k3z);

k4z=fz(x(j)+h,y(j)+h*k3y,z(j)+h*k3z);

y(j+h)=y(j)+h/6*(k1y+2*k2y+2*k3y+k4y);

z(j+h)=z(j)+h/6*(k1z+2*k2z+2*k3z+k4z);

end

x = 1:N:xfinal;

figure;

plot(x,y,z,'LineWidth',2);

xlabel('X');

ylabel('Y');

I took out the μ value and decreased the range to (0,10) just until I can understand how to form the code and this is what I have done but there is a problem with the plot. Can I get some help with it. Thanks.

Answer 1

[Torsten]

fy=@(x,y,z) z;
fz=@(x,y,z) (1-y^2)*z-y;
x(1)=0;
z(1)=0;
y(1)=1;
h=0.1;
xfinal=10;
N=ceil((xfinal-x(1))/h);
for j=1:N 
    x(j+1)=x(j)+h;
    k1y=fy(x(j),y(j),z(j));
    k1z=fz(x(j),y(j),z(j));
    k2y=fy(x(j)+h/2,y(j)+h/2*k1y,z(j)+h/2*k1z);
    k2z=fz(x(j)+h/2,y(j)+h/2*k1y,z(j)+h/2*k1z);
    k3y=fy(x(j)+h/2,y(j)+h/2*k2y,z(j)+h/2*k2z);
    k3z=fz(x(j)+h/2,y(j)+h/2*k2y,z(j)+h/2*k2z);
    k4y=fy(x(j)+h,y(j)+h*k3y,z(j)+h*k3z);
    k4z=fz(x(j)+h,y(j)+h*k3y,z(j)+h*k3z);
    y(j+1)=y(j)+h/6*(k1y+2*k2y+2*k3y+k4y);
    z(j+1)=z(j)+h/6*(k1z+2*k2z+2*k3z+k4z);
end
figure;
plot(x,[y;z],'LineWidth',2);
xlabel('X');
ylabel('Y');

[Cubii4]

Thank you so much.

Answer 2

[Cris LaPierre]

https://www.mathworks.com/matlabcentral/answers/?search_origin=ans_leaf&term=runge-kutta+4th+order

[Cubii4]

I still have some problems with the plot. I wrote the code in the comment above if I can get some help with it. Thanks.