function [x,y] = Runge_Kutta(x0,y0,b,n,f)
% Runge Kutta method with n subintervals to solve the initial value problem dy/dx=f(x,y) with
% initial conditions y(x0)=y0 on the interval [x0,b]. The function
% f(x,y) must be entered as a string. 

% For example, if you wish to solve dy/dx=y*cos(x) with y(0)=1 on the
% interval [0,15] with 50 intervals you would enter
% Runge_Kutta(0,1,15,50,'y.*cos(x)')

f = str2func(['@(x,y) ',f]); % converts the string input into a function handle

h=(b-x0)/n;
x=zeros(n+1,1); y= zeros(n+1,1); % Pre-allocoting vectors for speed
x(1)=x0; y(1)=y0;
for i=1:n
    x(i+1)=x0+i*h;
    k1=f(x(i),y(i));
    k2=f(x(i)+h/2,y(i)+h/2*k1);
    k3=f(x(i)+h/2,y(i)+h/2*k2);
    k4=f(x(i+1),y(i)+h*k3);
    y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6;
end
end

