Fortran programming to solve linear equation for ode

[ra_forever8]

Fint the exact solution of the system
dy/dt = -15y-25z
dz/dt=-47y-85z
with inital condition y(0)=2, z(0)=5
either by writing the equation in matrix form as dx/dt =AX where x=(y z) and diagonalising the matrix A, or otherwise.
Using fortran programming with second order adam bashforth method, plot the numerical soltuions for y_i and z_i against the exact solution between t=0 and t=0.1 for different time steps h(large and small), chosen to distinguish between instability and stability(or in the case of an unconditionally stable scheme, between being oscillatory and non oscillatory)
Here is my code

!Second order Adams-Bashforth method for ODE
!dy/dt= -15y-25z
!dz/dt=-47y-85z
! with intial condition y(0)=2, z(0)=5

Program adamstwo
Implicit None
Real, allocatable :: y(:),t(:), z(:), texact(:),yexact(:),zexact(:)
Real:: tend, h, k1,k2,k3,k4
Real, parameter :: exp=2.718282
Integer:: NI, i
Print*, 'Enter the final time '
read*, Tend
Print*, 'Enter number of timesteps to take'
read*, NI
h= Tend/NI
Print*, 'This gives stepsize h=',h
allocate (t(0:NI), y(0:NI),z(0:NI))
!Initial Conditions
t(0) = 0
y(0) = 2
z(0) = 5

!After using runge kutta method, I found out k1 =-155 and k2= 3860*h-155 ,
k1 = -155
k2 = 3860*h-155
!we know that y(n+1) =y(n) + h/2(k1+k2)at n=0
y(1) = y(0) + h/2 *( k1 + k2)

!After using runge kutta method, I found out k3 =-155 and k4= 3860*h-155
k3= -519
k4= 44068*h-519
!we know that z(n+1) =z(n) + h/2(k3+k4) at n=0
z(1) = z(0) + h/2 *( k3 + k4)

t(1) = h

open(10,file='adams_two results.m')
! Loop through the number of steps to calculate the following at each step
do i = 2, NI
t(i) = i*h

!Second order Adam bashforth for all n
y(i+1) = y(i)+ (h/2)*(3*(-15*y(i)-25*z(i))-(-15*y*(i-1)-25*z*(i-1)))
z(i+1)= z(i)+ (h/2)*(3*(-47*y(i)-85*z(i))-(-47*y*(i-1)-85*z*(i-1)))

end do

!Print out the Approximate solution
write(10,*) 'ApproximateSolution =[', t(0),y(0),z(0)
do i = 0, NI
write(10,*) t(i),y(i),z(i)
end do
write(10,*) t(NI), y(NI),z(NI),']'
allocate (texact(0:NI), yexact(0:NI),zexact(0:NI))
texact(0)=0
yexact(0)=2
zexact(0)=2
do i = 1, NI
texact(i) = i*h
yexact(i)=1/(4*(i*h)**4 +1)
yexact(i)= 0.4385* exp**((-50+20*sqrt(6.0))*(i*h)) + 1.5613 *exp**((-50-20*sqrt(6.0))*(i*h))
zexact(i)= -0.2455* exp**((-50+20*sqrt(6.0))*(i*h)) + 5.2453* exp**((-50-20*sqrt(6.0))*(i*h))
end do

!Print out the exact solution
write(10,*) 'ExactSolution = [',texact(0), yexact(0),zexact(0)
do i = 0, NI
write(10,*) texact(i), yexact(i),zexact(i)
end do
write(10,*) texact(NI), yexact(NI),zexact(NI),']'
write(10,*) "plot(ApproximateSolution(:,1),ApproximateSolution(:,2),'g',ExactSolution(:,1),ExactSolution(:,2),'r')"
write(10,*) "xlabel('time'),ylabel('y'),legend('Approximate AB[2] Solution','Exact Solution')"
close(10)

end program adamstwo

I got the run time error :
Error 112,Reference to undefined variable,array element or function result(/UNDEF)
main - in file adamstwo.f95 at line 44 [+0877] ---
(i don't have MATLAB file because of this runtime in fortran programming.)

Please help.
(Note:I have calculated the exact solution of the system by hand which is correct)

 

[Mark44]

Here are lines 44 and 45 of your code. I have included both lines, because even if you fix line 44, the compiler will still give you an error for line 45

y(i+1) = y(i)+ (h/2)*(3*(-15*y(i)-25*z(i))-(-15*y*(i-1)[/color]-25*z*(i-1)[/color])) 
 z(i+1)= z(i)+ (h/2)*(3*(-47*y(i)-85*z(i))-(-47*y*(i-1)[/color]-85*z*(i-1)[/color]))

There are no y and z variables in your code - there are arrays named y and z. Where you have y * <something> and z * <something> in the lines above, you should have y(i - 1) and z(i - 1), to access the element in each array at index i - 1.

It's possible that there are other errors, but if the fix the above two lines, that should take care of the run-time error you're seeing.

 

[ra_forever8]

ok i changed it
y(i+1) = y(i)+ (h/2)*(3*(-15*y(i)-25*z(i))-(-15*y(i-1)-25*z(i-1)))
z(i+1)= z(i)+ (h/2)*(3*(-47*y(i)-85*z(i))-(-47*y(i-1)-85*z(i-1)))
and added this code : t(NI)=NI *h and its working. Thank you mark 44 and others too.

 

[ra_forever8]

i want to plot two graph : yn against exact solution and zn against exact solution separately
Here is fortran code to open file for matlab:
write(10,*) "plot(yn(:,1),yn(:,2),'b',ExactSolution(:,1),ExactSolution(:,2),'r')"
write(10,*) "xlabel('time'),ylabel('y'),legend('yn','Exact Solution')"
write(10,*) "plot(zn(:,1),zn(:,2),'g',ExactSolution(:,1),ExactSolution(:,2),'r')"
write(10,*) "xlabel('time'),ylabel('y'),legend('zn','Exact Solution')"
But it does not work when i open the MATLAB file , it gives the numerical solution of yn, zn and exact solution. but when i run to plot the two graphs , it only plot one graph which is zn aginst exact solution. it does not open or plot yn against exact solution.

 

[ra_forever8]

!Second order Adams-Bashforth method for ODE
!dy/dt= -15y-25z
!dz/dt=-47y-85z
! with intial condition y(0)=2, z(0)=5
!Rabindra Gurung Qs part2
Program adamstwo
Implicit None
Real, allocatable :: y(:),t(:), z(:), texact(:),yexact(:),zexact(:)
Real:: tend, h, k1,k2,k3,k4
Real, parameter :: exp=2.718282
Integer:: NI, i
Print*, 'Enter the final time '
read*, Tend
Print*, 'Enter number of timesteps to take'
read*, NI
h= Tend/NI
Print*, 'This gives stepsize h=',h
allocate (t(0:NI), y(0:NI),z(0:NI))
!Initial Conditions
t(0) = 0.0
y(0) = 2.0
z(0) = 5.0

!After using runge kutta method, I found out k1 =-155 and k2= 3860*h-155 ,
k1 = -155.0
k2 = 3860.0*h-155.0
!we know that y(n+1) =y(n) + h/2(k1+k2)at n=0
y(1) = y(0) + h/2 *( k1 + k2)

!After using runge kutta method, I found out k3 =-155 and k4= 3860*h-155
k3= -519.0
k4= 44068.0*h-519.0
!we know that z(n+1) =z(n) + h/2(k3+k4) at n=0
z(1) = z(0) + h/2.0 *( k3 + k4)

open(10,file='adams_two results.m')
! Loop through the number of steps to calculate the following at each step
do i = 1, NI-1
t(i) = i*h

!Second order Adam bashforth for all n
y(i+1) = y(i)+ (h/2.0)*(3.0*(-15.0*y(i)-25.0*z(i))-(-15.0*y(i-1)-25.0*z(i-1)))
z(i+1)= z(i)+ (h/2.0)*(3.0*(-47.0*y(i)-85.0*z(i))-(-47.0*y(i-1)-85.0*z(i-1)))
end do
t(NI)=NI *h

!Print out the yn
write(10,*) 'yn =[', t(0),y(0)
do i = 0, NI
write(10,*) t(i),y(i)
end do
write(10,*) t(NI), y(NI),']'

!Print out the zn
write(10,*) 'zn =[', t(0),z(0)
do i = 0, NI
write(10,*) t(i),z(i)
end do
write(10,*) t(NI), z(NI),']'

allocate (texact(0:NI), yexact(0:NI),zexact(0:NI))
texact(0)=0
yexact(0)=2
zexact(0)=2
do i = 1, NI
texact(i) = i*h
yexact(i)= 0.4385* exp**((-50.0+20.0*sqrt(6.0))*(i*h)) + 1.5613 *exp**((-50.0-20.0*sqrt(6.0))*(i*h))
zexact(i)= -0.2455* exp**((-50.0+20.0*sqrt(6.0))*(i*h)) + 5.2453*exp**((-50.0-20.0*sqrt(6.0))*(i*h))
end do
!Print out the exact solution
write(10,*) 'ExactSolution = [',texact(0), yexact(0),zexact(0)
do i = 0, NI
write(10,*) texact(i), yexact(i),zexact(i)
end do
write(10,*) texact(NI), yexact(NI),zexact(NI),']'
write(10,*) "plot(yn(:,1),yn(:,2),'b',ExactSolution(:,1),ExactSolution(:,2),'r')"
write(10,*) "xlabel('time'),ylabel('y'),legend('yn','Exact Solution')"
write(10,*) "plot(zn(:,1),zn(:,2),'g',ExactSolution(:,1),ExactSolution(:,2),'r')"
write(10,*) "xlabel('time'),ylabel('y'),legend('zn','Exact Solution')"
write(10,*)"hold all"
close(10)
end program adamstwo

(Here is my final corrected fortran code and I would like to plot yn against exact solution and zn against exact solution separately but it giving only one plot in matlab.)