RK4 ode45 solver segmentation fault

[millertime1234]

Hey guys so I've been writing a small ode45 solver over the past couple of days and finally got it to compile and run. However after all of my values are read in by the user I get a segmentation fault. I'm pretty sure everything is allocated for correctly and I've been messing with this for a couple hours now. Any ideas? here's a copy of my code.


!ODE SOLVER
!number of steps should be large
!value of delx should be small
program odesolver
real(kind=8) :: x0, delx, x
real(kind=8), allocatable :: ysol,, xsol), k1), k2), k3), k4), y0), y)
integer :: N,Nstep,i,j
write(*,*),"enter the number of equations"
read(*,*), N
write(*,*),"enter the number of steps"
read(*,*), Nstep
write(*,*),"enter the value of delx"
read(*,*), delx
allocate(ysol(N,0:Nstep),xsol(0:Nstep),k1(N),k2(N),k3(N),k4(N),y0(N),y(N))
write(*,*),"enter the value of x0"
read(*,*), xo
DO i=1,N
write(*,*),"enter the values of Y0",i
read(*,*), YO
end DO

xsol(0)=x0
DO i=1,N
ysol(i,0)=y0(i)
end DO
x=x0
DO i=1,N
y(i)=ysol(i,0)
end DO
DO i=1,Nstep
call RK4(Nstep,N,delx,x,y,xsol,ysol,k1,k2,k3,k4)
xsol(i)=x
DO j=1,N
ysol(j,i)=y(j)
end DO
end DO
end program
subroutine RK4(Nstep,N,delx,x,y,xsol,ysol,k1,k2,k3,k4)
integer :: N,Nstep,i,j
real(kind=8) :: y(N), k1(N), k2(N), k3(N), k4(N), xsol(0:Nstep), ysol(N,0:Nstep), xshift, yshift(N)
call Derivs(x,y,k1)
xshift = x+.05*delx
DO j=1,N
yshift(j)=y(j)+0.5*k1(j)*delx
end DO
call Derivs(N,xshift,yshift,k2)
xshift=x+0.5*delx
DO j=1,N
yshift(j)=y(j)+0.5*k2(j)*delx
end DO
call Derivs(N,xshift,yshift,k3)
xshift = x+delx
DO j=1,N
yshift(j)=y(j)+k3(j)*delx
end DO
call Derivs(N,xshift,yshift,k4)
xshift = x+delx
DO j=1,N
y(j)=y(j)+(delx/6.0)*(k1(j)+2.0*k2(j)+2.0*k3(j)+k4(j))
end DO
end subroutine

subroutine Derivs(N,x,y,f)
real(kind=8) :: f(N)
f(1)=y(2)
f(2)=5.0*sin(x)-20*abs(y(2))*y(2)-8.0*y(1)*y(1)*y(1)
end subroutine

Any help would be appreciated!

 

[xwb]

Which compiler/OS are you using?

If Linux, build with -g -ggdb, then run the program with gdb. It should tell you where the fault is occuring.
If MS/Intel, run it on visual studio and it should tell you where it is getting the segv.

 

[salgerman]

Simple visual inspection reveals a few things wrong with the source.

The entering of values for Y0 is being done wrong...first, you wrote YO(<-letter 'O') and not Y0(<-number zero); second, if you are going to do it that way, you need to provide index, too...read(*,*) Y0(i)

In your subroutine RK4, the first call to sub Derivs only passes 3 arguments...it is missing the very first one.

In your subroutine Derivs, you did not declare the x and y parameters and while this may not be a problem with x, it is with y since it is an array!

 

[mikrom]

To ask the user for entering the number of the equations N is misleading, because it seems that you have strictly defined 2 equations in your program .
I suggest to use implicit none and declare all variables.