#####################################
# Integrate a first order ODE using
# 2nd order Runge-Kutta
#####################################

# Define the phase speed

f:= t-> cos(t):

##
# Set initial condition and step length and number of steps
##

y[0]:=1:  h:= 0.1:  N:= 100:

##
# Advance solution using N Runge-Kutta steps
##


t[0]:=0:
for i from 1 to N do

###
# Step from t[i-1] to t[i] using 
# the provisional step tt inbetween
###

# 1. Compute f(t,y)

	f0:= f(t[i-1]);

# 2. Set provisional value yt

	yt:= y[i-1] + h*f0;

# 3. Compute phase speed ft at provisional point

	ft:= f(t[i-1]+h);

# 4. Compute final phase speed ff

	ff:= (1/2)*f0 + (1/2)*ft;

# 5. Set new solution value at t[i]

	y[i]:= y[i-1] + h*ff;

# Update time

	t[i]:= t[i-1] + h;

end do;


##
# Build set of data points for plotting
##

l:= [[t[n],y[n]] $n=0..N]:

##
# Plot approximate solution against exact solution
##

with(plots):
plotsetup(x11):

#plot(l, style=line,symbol=circle);

F:= plot(l, style=point,symbol=circle):
G:= plot(1+sin(t),t=0..t[N],colour=blue):
display({F,G});