Runge Kutta 4ª ordem

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

clear all
clc
syms y
syms x
syms z
syms w
syms p1
syms p2
syms p3
p1 = input('Insira a equação diferencial dy/dx: ');
p2 = input('Insira a equação diferencial dz/dx: ');
p3 = input('Insira a equação diferencial dw/dx: ');
h = input('Insira o valor do passo: ');
a = input('Insira o valor de x inicial: ');
y0 = input('Insira o valor de y inicial: ');
z0 = input('Insira o valor de z inicial: ');
w0 = input('Insira o valor de w inicial: ');
b = input('Insira o valor de x final: ');
n = (b-a)/h;
i=0;
fprintf('\n\t a \t\t y \t\t\t z \t\t\t w \t\t\n')
for i=1:n
k1 =subs(p1,[x,y,z,w],[a,y0,z0,w0]);
o1 =subs(p2,[x,y,z,w],[a,y0,z0,w0]);
m1 =subs(p3,[x,y,z,w],[a,y0,z0,w0]);
k2 = subs(p1,[x,y,z,w],[a+0.5*h,y0+0.5*k1*h,z0+0.5*o1*h,w0+0.5*m1*h]);
o2 = subs(p2,[x,y,z,w],[a+0.5*h,y0+0.5*k1*h,z0+0.5*o1*h,w0+0.5*m1*h]);
m2 = subs(p3,[x,y,z,w],[a+0.5*h,y0+0.5*k1*h,z0+0.5*o1*h,w0+0.5*m1*h]);
k3 = subs(p1,[x,y,z,w],[a+0.5*h,y0+0.5*k2*h,z0+0.5*o2*h,w0+0.5*m2*h]);
o3 = subs(p2,[x,y,z,w],[a+0.5*h,y0+0.5*k2*h,z0+0.5*o2*h,w0+0.5*m2*h]);
m3 = subs(p3,[x,y,z,w],[a+0.5*h,y0+0.5*k2*h,z0+0.5*o2*h,w0+0.5*m2*h]);
k4 = subs(p1,[x,y,z,w],[a+h,y0+k3*h,z0+o3*h,w0+m3*h]);
o4 = subs(p2,[x,y,z,w],[a+h,y0+k3*h,z0+o3*h,w0+m3*h]);
m4 = subs(p3,[x,y,z,w],[a+h,y0+k3*h,z0+o3*h,w0+m3*h]);
a = a+h;
y0 = y0+1/6*(k1+2*k2+2*k3+k4)*h;
y0 = double(y0);
z0 = z0+1/6*(o1+2*o2+2*o3+o4)*h;
z0 = double(z0);
w0 = w0+1/6*(m1+2*m2+2*m3+m4)*h;
w0 = double(w0);
formatSpec = ' %.3f %.8f %.8f %.8f\n';
fprintf(formatSpec,a,y0,z0,w0)
end