Runge Kutta 44 en Matlab

clear all;close all;clc
fun=input('Ingresa la función f(t,y) entre comillas: ');
f=inline(fun,'t','y');
y(1)=input('Ingrese la condición inicial y(0): ');
t0=input('Ingrese t0: ');
max=input('Ingrese el t final: ');
h=input('Ingrese el paso: ');
t=t0:h:max;
n=length(t);


for i=1:n-1
    
    k1=feval(f,t(i),y(i));
    k2=feval(f,t(i)+h/2,y(i)+h*k1/2);
    k3=feval(f,t(i)+h/2,y(i)+h*k2/2);
    k4=feval(f,t(i)+h,y(i)+h*k3);
    
    y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);
    
end


plot(t,y)