function [rk4T, rk4W] = rk4(Fun, t0, tlast, w0, N)
%%%%%%%%%%% Runge-Kutta Method of Order 4 for Initial Value Problems %%%%%
%% 
%% Solve y'= f(t,y),   a<= t <=b,  y(a)=  alpha
%%
%% It is assumed that you have the function f(t,y) defined elsewhere 
%% as a function in a .m file
% INPUT:
% Fun   - String containing name of function m-file (enclose in '  ')
% t0    - Initial value of t.
% tlast - Last value of t.
% w0    - Initial value.
% N     - The number of steps
%
% OUTPUT:
% rk4T  - time points (column-vector).
% rk4W  - Solution column-vector per rk4T-values.
%
% The result can be displayed by: plot(rk4T, rk4W).
%% Step 1
h=(tlast-t0)/N;t=t0;w=w0;
rk4T=t;rk4W=w;

%% Step 2
for i=1:N;
 k1= h*feval(Fun,t,w);  
 k2= h*feval(Fun,t+(h/2),w+(k1/2));  
 k3= h*feval(Fun,t+(h/2),w+(k2/2));  
 k4= h*feval(Fun,t+h,w+k3);  
 w = w + (k1+2*k2+2*k3+k4)/6; %%Step 3
 t = t0 + i*h;                      %%Step 3
 rk4T=[rk4T;t]; rk4W = [rk4W;w];               %%Step 4
end; 

