Runge-Kutta ordem 4

Boa tarde Galera!!!

bom quero neste post atualizar o metodo de runge kutta, pois o primeiro que postei estava faltando algumas linhas de comando e outras linhas estavam incorretas... entao quero agradecer aos meus amigos da eng eletrica q ajudaram a corrigir o codigo!!!
ate onde testei o metodo funciona para qualquer funcao.

syms x y
f  = input ('Entre com a equação: ');
x0 = input ('Entre com o valor inicial para x: ');
xf = input ('entre com o valor final para x');
y0 = input ('Entre com o valor inicial y: ');
h  = input ('Entre com o passo desejado: ');
n  = (xf - x0)/h;

for i= 1:n
    k1 = h * subs(f,{x,y}, {x0,y0});
    k2 = h * subs(f,{x,y},{x0 + h/2, y0 +k1/2});
    k3 = h * subs(f,{x,y},{x0 + h/2, y0 +k2/2});
    k4 = h * subs(f,{x,y},{x0 + h, y0 +k3});
    ya= y0+ (1/6) * (k1+(2*k2)+(2*k3)+k4);
    x0 = x0 + h;
    y0 = ya;
    fprintf('y%.0f= %.10f.\n',n,ya)
   
end

Entao esta ai... agora certinhu... este projeto valeu 10 pro prof de cal numerico, e espero que ajude vcs...

beijinhus
adoro vcs