un: matrix for storing solution at every time step
for j = 2:nt+1
rhs(nx,dx,un,r,alpha)
for i = 2:nx
ut[i] = un[i] + dt*r[i] # 1st step
end
rhs(nx,dx,ut,r,alpha)
for i = 2:nx
ut[i] = 0.75*un[i] + 0.25*ut[i] + 0.25*dt*r[i] # 2nd step
end
rhs(nx,dx,ut,r,alpha)
for i = 2:nx
un[i] = (1.0/3.0)*un[i] + (2.0/3.0)*ut[i] + (2.0/3.0)*dt*r[i] # 3rd step
end
k = k+1 # index for solution storage at every time step
u[:,k] = un[:]
end