# Help with Runge-Kutta 4

**URL:** <https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910>\
**Category:** Homework\
**Created:** [December 3, 2023, 4:15pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910 "2023-12-03T16:15:08Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![4ulalia](https://avatars.discourse-cdn.com/v4/letter/4/df705f/32.png) [@4ulalia](https://fortran-lang.discourse.group/u/4ulalia)\
**Post date:** [December 3, 2023, 4:15pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/1 "2023-12-03T16:15:08Z")

</div>

Hello, I’m trying to program a subroutine that solves differential equations using the Runge-Kutta 4 method. To try if it’s correct I’m implementing it on a pendulum equation so I can plot it and check if the result is correct. Although i though it was a very simple formula, I am having some struggles because when excecuting I have the following message : Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

This messages somewtimes occur and I tried to search for the memory error but can seem to find it on the code. Can someone give it a go?

Thanks!

`````auto
PROGRAM PREPA8
	IMPLICIT NONE
	DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud
	INTEGER:: N, i, iter, nequ
	DOUBLE PRECISION:: yin_0(2), yy1(2)
	EXTERNAL:: EDO, pendolSimple
	COMMON/DADES/V

	!Dades ctts::
	V=-2.4d0
	pi=4d0*atan(1d0)
	grav = 3.71d0 !m/s^2
	longitud = 0.45d0 !m

	n=1300 !número de passos
	!Inicialitzem els valors amb les condicions inicials
	yin_0(1) = 0.02d0
	yin_0(2) = 0d0
	dx=0.0003d0

	do i = 1, N
		call RungeKutta4order(nequ, dx, yin_0,yy1, pendolSimple)
		write(*,*) dx*i, yy1
		!canviem les variables pel bucle
		yin_0=yy1
	enddo
	
END PROGRAM PREPA8

	IMPLICIT NONE
	INTEGER:: nequ, i
	DOUBLE PRECISION:: yyin(nequ), yyout(nequ), k1(nequ), k2(nequ), k3(nequ), k4(nequ), yaux(nequ)
	DOUBLE PRECISION:: x, dx

	x=0d0 !realment no necessitem x

	call func(nequ, x, yyin, k1) !k1 = f(x0, y0)
	yaux=yyin+(dx+k1/2d0)
	call func(nequ, x+(dx/2d0), yaux, k2)
	yaux=yyin +(dx*k2/2d0)
	call func(nequ, x+(dx/2d0), yaux, k3)
	yaux=yyin+(dx*k3)
	call func(nequ, x+dx, yaux, k4)

	yyout=yyin+(dx/6d0)*(k1+k2*2d0+2d0*k3+k4)

END SUBROUTINE RungeKutta4order

SUBROUTINE pendolSimple(yin, f)
	IMPLICIT NONE
	DOUBLE PRECISION:: yin(2), f(2), grav, longitud, massa
	grav = 3.71d0 !m/s^2
	longitud = 0.45d0 !m
	f(1) = yin(2)
	f(2) = -(grav/longitud) * sin(yin(1))
	return

END SUBROUTINE pendolSimple ````
`````

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [December 3, 2023, 5:56pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/2 "2023-12-03T17:56:20Z")

</div>

Hard to say since this example isn’t complete.

- you aren’t defining `nequ` anywhere
- set some compiler flags to warn you of uninitialized variables.
- put all your routines in a module
- use [a runge-kutta library](https://github.com/jacobwilliams/rklib) instead of reinventing the wheel

---

<div class="post-metadata">

**Author:** ![4ulalia](https://avatars.discourse-cdn.com/v4/letter/4/df705f/32.png) [@4ulalia](https://fortran-lang.discourse.group/u/4ulalia)\
**Post date:** [December 3, 2023, 6:19pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/3 "2023-12-03T18:19:14Z")

</div>

What do you mean it isn’t complete?  
I never use modules and they specifically asked me to do a Runge Kutta subroutine. I inicialized `nequ` but still getting the error

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:  
#0 0x101e3aafe  
#1 0x101e39cdd  
#2 0x7fff2041dd7c  
#3 0x101e28659  
#4 0x101e28978  
#5 0x101e28adf  
zsh: segmentation fault ./a.out

```PROGRAM
	IMPLICIT NONE
	DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud
	INTEGER:: N, i, iter
	INTEGER, PARAMETER:: nequ=2
	DOUBLE PRECISION:: yin_0(nequ), yy1(nequ),
	EXTERNAL:: pendolSimple
	COMMON/DADES/V

	!Dades ctts::
	V=-2.4d0
	pi=4d0*atan(1d0)
	grav = 3.71d0 !m/s^2
	longitud = 0.45d0 !m

	n=1300 !número de passos
	!Inicialitzem els valors amb les condicions inicials
	yin_0(1) = 0.02d0
	yin_0(2) = 0d0
	dx=0.0003d0

	do i = 1, N
		call RungeKutta4order(2, dx, yin_0,yy1, pendolSimple)
		write(*,*) dx*i, yy1
		!canviem les variables pel bucle
		yin_0 = yy1
	enddo
	
END PROGRAM PREPA8

SUBROUTINE RungeKutta4order(nequ, dx,yyin,yyout, func)
	IMPLICIT NONE
	INTEGER:: nequ, i
	DOUBLE PRECISION:: yyin(nequ), yyout(nequ), k1(nequ), k2(nequ), k3(nequ), k4(nequ), yaux(nequ)
	DOUBLE PRECISION:: x, dx

	x=0d0 !realment no necessitem x

	call func(nequ, x, yyin, k1) !k1 = f(x0, y0)
	yaux=yyin+(dx+k1/2d0)
	call func(nequ, x+(dx/2d0), yaux, k2)
	yaux=yyin +(dx*k2/2d0)
	call func(nequ, x+(dx/2d0), yaux, k3)
	yaux=yyin+(dx*k3)
	call func(nequ, x+dx, yaux, k4)

	yyout=yyin+(dx/6d0)*(k1+k2*2d0+2d0*k3+k4)

END SUBROUTINE RungeKutta4order

SUBROUTINE pendolSimple(yin, f)
	IMPLICIT NONE
	DOUBLE PRECISION:: yin(2), f(2), grav, longitud, massa
	grav = 3.71d0 !m/s^2
	longitud = 0.45d0 !m
	f(1) = yin(2)
	f(2) = -(grav/longitud) * sin(yin(1))
	return

END SUBROUTINE pendolSimple

```

---

<div class="post-metadata">

**Author:** ![NdotA](https://avatars.discourse-cdn.com/v4/letter/n/f14d63/32.png) [@NdotA](https://fortran-lang.discourse.group/u/NdotA)\
**Post date:** [December 3, 2023, 6:58pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/4 "2023-12-03T18:58:43Z")

</div>

What strikes me with your code:

Your Subroiutine pendolSimple is defined with two paramters but you call it from your RungeKutta subroutine with four.

Ah, and by the way, if grav is the Earth’s gravity you should check this value. And I doubt you would need double precision for this task, especially when you run it on a 64 bit machine.  
Cheers  
Norbert

---

<div class="post-metadata">

**Author:** ![4ulalia](https://avatars.discourse-cdn.com/v4/letter/4/df705f/32.png) [@4ulalia](https://fortran-lang.discourse.group/u/4ulalia)\
**Post date:** [December 3, 2023, 7:05pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/5 "2023-12-03T19:05:49Z")

</div>

Thank you, it was indeed the fact the number of parameters.

---

<div class="post-metadata">

**Author:** ![aledinola](https://avatars.discourse-cdn.com/v4/letter/a/3ec8ea/32.png) [@aledinola](https://fortran-lang.discourse.group/u/aledinola)\
**Post date:** [December 3, 2023, 8:03pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/6 "2023-12-03T20:03:15Z")

</div>

I think an advantage of putting all functions and subroutines into a module is that they will have an explicit interface and the compiler will check that you are calling the function with the correct arguments. What I found weird of your example is that you got this mistake only at run time and the compiler was not able to flag this during compilation.  
Using modules makes your code easier to debug and to find mistakes

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [December 3, 2023, 8:12pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/7 "2023-12-03T20:12:49Z")

</div>

I think the code you are pasting here is not the code you are compiling. A compiler would have caught some of these syntax errors (such as the missing `subroutine RungeKutta4order` declaration in the first one, or the extra comma in the second one).

---

<div class="post-metadata">

**Author:** ![Pap](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pap/32/1305_2.png) [@Pap](https://fortran-lang.discourse.group/u/Pap)\
**Post date:** [December 3, 2023, 9:07pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/8 "2023-12-03T21:07:39Z")

</div>

Obviously, this is a student’s homework. Which is ok, but apparently you missed at least one line of code while copy/pasting.

> [@4ulalia](#):
>
> I never use modules and they specifically asked me to do a Runge Kutta subroutine.

The fact they asked you to do a subroutine doesn’t mean said subroutine can’t be in a module.  
Also, if you want to learn Fortran, you really need to learn how to use modules, so that old fashion statements (which you use in the code, like `common` and `external`) won’t be needed. The code will be much easier to understand (and expand later on.)

Fortran is much easier to learn compared to other languages like C or C++, but you still need to put some effort into learning it…

---

<div class="post-metadata">

**Author:** ![4ulalia](https://avatars.discourse-cdn.com/v4/letter/4/df705f/32.png) [@4ulalia](https://fortran-lang.discourse.group/u/4ulalia)\
**Post date:** [December 4, 2023, 10:19am UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/9 "2023-12-04T10:19:59Z")

</div>

Thank you for the advise, will start to use modules from now on.  
I have encountered another problem while creating a subroutine that is supposed to return the excited level of a particle. To do so I created a subroutine `integrate` that integrates the Schrodinger equation using Runge-Kutta4 but it seems that it’s not updating how it’s suppose to be yyin and yyout at each iteration. Does anyone know what the problem might be?

```auto
SUBROUTINE integrate(nequ, E0, a, b, N, yyin, phi)
	IMPLICIT NONE
	DOUBLE PRECISION:: E0, phi, a, b
	DOUBLE PRECISION:: dx, yyin(nequ), yyout(nequ), local_yyin(nequ)
	INTEGER:: nequ, i, N
	EXTERNAL:: EDO

	dx=(b-a)/N
	local_yyin=yyin
	DO i = 1, N
		call RungeKutta4order(nequ, dx, local_yyin, yyout, edo)
		print*, dx*i, yyin, yyout
		local_yyin=yyout
	ENDDO

	phi=local_yyin(1)

END SUBROUTINE integrate

SUBROUTINE secant(E1, E2, phi1, phi2, E3)
	IMPLICIT NONE
	DOUBLE PRECISION:: E1, E2, phi1, phi2, E3

	E3=(E1 * phi2 - E2 * phi1) / (phi2 - phi1)
END SUBROUTINE secant

SUBROUTINE tir(nequ, E1, E2, error, a,b,N, yyin)
	IMPLICIT NONE
	DOUBLE PRECISION:: E1, E2, E3, a, b
	DOUBLE PRECISION:: phi1, phi2, phi3, dx, x, error, yyin(nequ), yyout(nequ)
	INTEGER:: i, nequ, N, iter

	i=0 !incialitzo
	iter=500
	!Pas 1 i 2

	do i = 1, iter
		call integrate(nequ, E1, a, b, N,yyin, phi1)
		call integrate(nequ, E2, a, b, N, yyin, phi2)
		print*, phi1, phi2

		call secant(E1, E2, phi1, phi2, E3)
		call integrate(nequ, E3, a, b, N,yyin, phi3)

		if (abs(phi3).LT.error) then
			write(*,*) 'Ha convergit', i, E3
			exit

	    else
	    	E1 = E2
		    E2 = E3
		end if
	end do	

END SUBROUTINE tir
SUBROUTINE RungeKutta4order(nequ, dx, yyin, yyout, func)
	IMPLICIT NONE
	INTEGER:: nequ, i
	DOUBLE PRECISION, DIMENSION(nequ):: yyin, yyout, k1, k2, k3, k4, yaux
	DOUBLE PRECISION:: x, dx

	call func(nequ, x, yyin, k1) !k1 = f(x0, y0)

	yaux=yyin+(dx*k1/2d0)
	x=x+0.5d0*dx
	call func(nequ, x, yaux, k2) !k2
	
	yaux=yyin +(dx*k2/2d0)
	x=x+0.5d0*dx
	call func(nequ, x, yaux, k3) !k3
	
	yaux=yyin+(dx*k3)
	x=x+0.5d0*dx
	call func(nequ, x+dx, yaux, k4) !k4

	yyout=yyin+(dx/6d0)*(k1+k2*2d0+2d0*k3+k4)
	
END SUBROUTINE RungeKutta4order

!L'edos en aquest cas és equivalent a la del pèndol simple en aproximacions petites.
!yin=(phi, dphi/dx)
!dyout=(dphi/dx, d^2phi/dx^2)
SUBROUTINE EDO(nequ, x, yin, dyout)
	IMPLICIT NONE
	DOUBLE PRECISION:: x, yin(nequ), dyout(nequ), E, V
	INTEGER:: nequ
	COMMON/ENERGIES/E, V

	dyout(1)= yin(2)
	dyout(2) = 2D0*(V-E)* yin(1) !Correspon a la segona derivada que ens dona l'equació de Schrodinger

END SUBROUTINE EDO

PROGRAM PREPA8
	IMPLICIT NONE
	DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud, a, b
	INTEGER:: N, i, iter
	INTEGER, PARAMETER:: nequ=2
	DOUBLE PRECISION:: yyin(nequ), yyout(nequ), phi(200), phi_f(2)
	EXTERNAL:: EDO, pendolSimple
	COMMON/DADES/V, pi

	OPEN(10, file='prepa8.dat')
	!Dades ctts::
	V=-2.4d0
	pi=4d0*atan(1d0)
	
	yyin(1)=0.02d0
	yyin(2)=0.0d0
	E1 = 1d0
	E2= 3d0
	error=10-5d0
	call tir(nequ, E1, E2, error, 0d0,1d0,500, yyin)
END PROGRAM PREPA8 

```

---

<div class="post-metadata">

**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\
**Post date:** [December 4, 2023, 10:35am UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/10 "2023-12-04T10:35:58Z")

</div>

> [@NdotA](#):
>
> And I doubt you would need double precision for this task, especially when you run it on a 64 bit machine.

The use of `double precision` or not does not depend on the architecture of the machine.

---

<div class="post-metadata">

**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\
**Post date:** [December 4, 2023, 10:48am UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/11 "2023-12-04T10:48:41Z")

</div>

> [@4ulalia](#):
>
> Does anyone know what the problem might be?

At least, the variables in the `ENERGY` common block are not assigned any value.

Please, do yourself a favor: stop using common blocks and put everything in a module.

---

<div class="post-metadata">

**Author:** ![NdotA](https://avatars.discourse-cdn.com/v4/letter/n/f14d63/32.png) [@NdotA](https://fortran-lang.discourse.group/u/NdotA)\
**Post date:** [December 4, 2023, 11:26am UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/12 "2023-12-04T11:26:33Z")

</div>

Anyway, even 32 bit precision should be sufficient to check the outcome. But of course ‘DOUBLE PRECISION’ does not do any harm.

---

<div class="post-metadata">

**Author:** ![Pap](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pap/32/1305_2.png) [@Pap](https://fortran-lang.discourse.group/u/Pap)\
**Post date:** [December 4, 2023, 12:03pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/13 "2023-12-04T12:03:31Z")

</div>

Normally double precision isn’t needed even for a double pendulum with a friction factor added in the equations - unless someone needs insane precision, I guess.

Note that in this case we have a fixed step Runge-Kutta method - which is fine for educational purposes, but should never be used for real-world problems. Obviously this is student’s homework for some class though, and a fixed step Runge-Kutta is used to introduce the general idea of such methods.

---

<div class="post-metadata">

**Author:** ![4ulalia](https://avatars.discourse-cdn.com/v4/letter/4/df705f/32.png) [@4ulalia](https://fortran-lang.discourse.group/u/4ulalia)\
**Post date:** [December 4, 2023, 2:48pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/14 "2023-12-04T14:48:03Z")

</div>

I do define the variables in the COMMON at my main program (at the end of the code). I was told to use COMMONS but it’s definitely not necessary in this case so I will delete it. I still have the problem at the subroutine integrate where the values I print during the loop are always the initial values, they don’t get updated even though I wrote `yyin=yyout` at the end of the loop. I believe this is why the subroutine `tir` doesn’t work.  
I am learning to program and debugging my own codes but it can get very frustrating so help and advice is always great.  
Can someone see the problem with the code?? I really would appreciate some help

---

<div class="post-metadata">

**Author:** ![Pap](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pap/32/1305_2.png) [@Pap](https://fortran-lang.discourse.group/u/Pap)\
**Post date:** [December 4, 2023, 5:04pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/15 "2023-12-04T17:04:01Z")

</div>

The first thing you should do, and do it religiously in any program, is to _always_ add the `intent` attribute in _all_ subroutine/function arguments see, e.g., [this](https://fortran-lang.org/en/learn/quickstart/organising_code/) for a quick description of the attribute.)

For example, in your integrate subroutine, I guess `phi` is `intent(out)` - or is it `intent(inout)`? Do yourself a favor and learn to use `intent`. Not only it makes your code more readable, but doing that alone may help you find where the problem is. Do it, and I’ll have a better look in the code.

Also, when you ask for help, please try to post a minimal program that reproduces the problem. Posting long code makes people way less willing to help.

---

<div class="post-metadata">

**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\
**Post date:** [December 4, 2023, 6:51pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/16 "2023-12-04T18:51:43Z")

</div>

> [@4ulalia](#):
>
> I was told to use COMMONS

Whoever told you that in 2023 has a few decade to catch up with Fortran…

> [@4ulalia](#):
>
> I still have the problem at the subroutine integrate where the values I print during the loop are always the initial values, they don’t get updated even though I wrote `yyin=yyout` at the end of the loop.

You are printing `yyin`, but your update is `yyin_local=yyout`. I guess you want to print `yyin_local`

---

<div class="post-metadata">

**Author:** ![4ulalia](https://avatars.discourse-cdn.com/v4/letter/4/df705f/32.png) [@4ulalia](https://fortran-lang.discourse.group/u/4ulalia)\
**Post date:** [December 5, 2023, 10:02am UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/17 "2023-12-05T10:02:58Z")

</div>

```auto
PROGRAM PREPA8
	IMPLICIT NONE
	DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud, a, b
	INTEGER:: N, i, iter
	INTEGER, PARAMETER:: nequ=2
	DOUBLE PRECISION:: yyin(nequ), yyout(nequ), phi(200), phi_f(2)
	EXTERNAL:: EDO, pendolSimple
	COMMON/DADES/V, pi

	OPEN(10, file='prepa8.dat')
	!Dades ctts::
	V=-2.4d0
	pi=4d0*atan(1d0)
	
	yyin(1)=0.02d0
	yyin(2)=0.0d0
	E1 = 1d0
	E2= 1.1d0
	error=10-5d0
	call tir(nequ, E1, E2, error, 0d0,1d0,500, yyin)

	CLOSE(10)

END PROGRAM PREPA8

! Subroutine for numerical integration using the Runge-Kutta 4th order method
! Purpose: Numerically integrates a system of ordinary differential equations (ODEs)
! Inputs:
! nequ: Number of equations in the system
! E0: Initial value for a parameter (input)
! a, b: Integration limits (input)
! N: Number of integration steps (input)
! yyin: Initial values of the dependent variables (input)
! Outputs:
! phi: Final value of the first dependent variable after integration (output)
SUBROUTINE integrate(nequ, E0, a, b, N, yyin, phi)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: nequ, N
  DOUBLE PRECISION, INTENT(IN) :: E0, a, b
  DOUBLE PRECISION, INTENT(IN) :: yyin(nequ)
  DOUBLE PRECISION, INTENT(OUT) :: phi
  DOUBLE PRECISION :: dx, yyout(nequ)
  INTEGER :: i
  EXTERNAL :: EDO, RungeKutta4order ! External function for ODE and Runge-Kutta method

  ! Calculate the step size for integration
  dx = (b - a) / N

  ! Perform numerical integration using Runge-Kutta 4th order method
  DO i = 1, N
    CALL RungeKutta4order(nequ, dx, yyin, yyout, edo)
    WRITE(*,*) dx * i, yyout
    yyin = yyout
  ENDDO

  ! Set the final value of the first dependent variable
  phi = yyin(1)
END SUBROUTINE integrate

! Subroutine for the secant method to find the root of a function
! Purpose: Implements the secant method
! Inputs:
! E1, E2: Initial guesses for the root (input)
! phi1, phi2: Function values corresponding to the initial guesses (input)
! Outputs:
! E3: Improved estimate for the root (output)
SUBROUTINE secant(E1, E2, phi1, phi2, E3)
  IMPLICIT NONE
  DOUBLE PRECISION, INTENT(IN) :: E1, E2, phi1, phi2
  DOUBLE PRECISION, INTENT(OUT) :: E3

  ! Apply the secant method formula to improve the root estimate
  E3 = (E1 * phi2 - E2 * phi1) / (phi2 - phi1)
END SUBROUTINE secant

! Subroutine for solving a nonlinear equation using the secant method and numerical integration
! Purpose: Solves a nonlinear equation iteratively
! Inputs:
! nequ: Number of equations in the system (input)
! E1, E2: Initial guesses for the root (input)
! error: Tolerance for convergence (input)
! a, b: Integration limits (input)
! N: Number of integration steps (input)
! yyin: Initial values of the dependent variables (input)
SUBROUTINE tir(nequ, E1, E2, error, a, b, N, yyin)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: nequ, N
  DOUBLE PRECISION, INTENT(IN) :: E1, E2, error, a, b
  DOUBLE PRECISION, INTENT(INOUT) :: yyin(nequ) ! Changed to INOUT for modification inside the subroutine
  INTEGER :: i, iter

  ! Initialize counters
  i = 0
  iter = 500

  ! Iteratively apply the secant method along with numerical integration
  DO i = 1, iter
    CALL integrate(nequ, E1, a, b, N, yyin, phi1)
    CALL integrate(nequ, E2, a, b, N, yyin, phi2)

    ! Display intermediate results
    PRINT*, 'Iteration:', i, 'E1:', E1, 'E2:', E2, 'phi1:', phi1, 'phi2:', phi2

    ! Apply the secant method to get an improved root estimate
    CALL secant(E1, E2, phi1, phi2, E3)
    CALL integrate(nequ, E3, a, b, N, yyin, phi3)

    ! Check for convergence
    IF (ABS(phi3) < error) THEN
      WRITE(*,*) 'Converged after', i, 'iterations. Root:', E3
      EXIT
    ELSE
      ! Update guesses for the root
      E1 = E2
      E2 = E3
    ENDIF
  ENDDO
END SUBROUTINE tir

```

These are the subroutines I’m using alongside the Runge-Kutta and Schrodinger diferential equation(EDO). I added a print at the `integrate` subroutine to check on the values. I’m getting all the time the same values (the intial ones) yyout even though I defined the update `yyin=yyout`. I think Runge-Kutta is okay because I tried it with a pendulum equation and it worked just fine. Does anyone see the problem

---

<div class="post-metadata">

**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\
**Post date:** [December 5, 2023, 11:13am UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/18 "2023-12-05T11:13:09Z")

</div>

> [@4ulalia](#):
>
> These are the subroutines I’m using alongside the Runge-Kutta and Schrodinger diferential equation(EDO). I added a print at the `integrate` subroutine to check on the values. I’m getting all the time the same values (the intial ones) yyout even though I defined the update `yyin=yyout`.

Sorry, but the code you have posted in your last message won’t compile at all, so it can’t be the code you are running. If you want some help try posting **the exact code** you are running.

---

<div class="post-metadata">

**Author:** ![aledinola](https://avatars.discourse-cdn.com/v4/letter/a/3ec8ea/32.png) [@aledinola](https://fortran-lang.discourse.group/u/aledinola)\
**Post date:** [December 5, 2023, 11:58am UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/19 "2023-12-05T11:58:08Z")

</div>

@4ulalia I would suggest that you put the main program and all functions and subroutines that you are using for this project in a single post here, adding also the instructions for compilation. Make sure that the code that you post here compiles and gives the same results. Then we can help 🙂

---

<div class="post-metadata">

**Author:** ![4ulalia](https://avatars.discourse-cdn.com/v4/letter/4/df705f/32.png) [@4ulalia](https://fortran-lang.discourse.group/u/4ulalia)\
**Post date:** [December 5, 2023, 2:13pm UTC](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910/20 "2023-12-05T14:13:19Z")

</div>

```auto

PROGRAM PREPA8
	IMPLICIT NONE
	DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud, a, b
	INTEGER:: N, i, iter
	INTEGER, PARAMETER:: nequ=2
	DOUBLE PRECISION:: yyin(nequ), yyout(nequ), phi(200), phi_f(2)
	EXTERNAL:: EDO, pendolSimple
	COMMON/DADES/V, pi

	OPEN(10, file='prepa8.dat')
	!Dades ctts::
	V=-2.4d0
	pi=4d0*atan(1d0)
	
	yyin(1)=0.02d0
	yyin(2)=0.0d0
	E1 = 1d0
	E2= 1.1d0
	error=10-5d0
	call tir(nequ, E1, E2, error, 0d0,1d0,500, yyin)

CLOSE(10)

END PROGRAM PREPA8

!SUBRUTINE RUNGE-KUTTA D'ORDRE 4
SUBROUTINE RungeKutta4order(nequ, dx, yyin, yyout, func)
	IMPLICIT NONE
	INTEGER:: nequ, i
	DOUBLE PRECISION, DIMENSION(nequ):: yyin, yyout, k1, k2, k3, k4, yaux
	DOUBLE PRECISION:: x, dx

	call func(nequ, x, yyin, k1) !k1 = f(x0, y0)

	yaux=yyin+(dx*k1/2d0)
	x=x+0.5d0*dx
	call func(nequ, x, yaux, k2) !k2
	
	yaux=yyin +(dx*k2/2d0)
	x=x+0.5d0*dx
	call func(nequ, x, yaux, k3) !k3
	
	yaux=yyin+(dx*k3)
	x=x+0.5d0*dx
	call func(nequ, x+dx, yaux, k4) !k4

	yyout=yyin+(dx/6d0)*(k1+k2*2d0+2d0*k3+k4)
	
END SUBROUTINE RungeKutta4order

!L'edos en aquest cas és equivalent a la del pèndol simple en aproximacions petites.
!yin=(phi, dphi/dx)
!dyout=(dphi/dx, d^2phi/dx^2)
SUBROUTINE EDO(nequ, x, yin, dyout)
	IMPLICIT NONE
	DOUBLE PRECISION:: x, yin(nequ), dyout(nequ), E, V
	INTEGER:: nequ
	COMMON/ENERGIES/E, V

	dyout(1)= yin(2)
	dyout(2) = 2D0*(V-E)* yin(1) !Correspon a la segona derivada que ens dona l'equació de Schrodinger

END SUBROUTINE EDO

UBROUTINE integrate(nequ, E0, a, b, N, yyin, phi)
	IMPLICIT NONE
	DOUBLE PRECISION:: E0, phi, a, b
	DOUBLE PRECISION:: dx, yyin(nequ), yyout(nequ)
	INTEGER:: nequ, i, N
	EXTERNAL:: EDO

	dx=(b-a)/N
	DO i = 1, N
		call RungeKutta4order(nequ, dx, yyin, yyout, edo)
		write(*,*) dx*i, yyout
		yyin=yyout
	ENDDO

	phi=yyin(1)

END SUBROUTINE integrate

SUBROUTINE secant(E1, E2, phi1, phi2, E3)
	IMPLICIT NONE
	DOUBLE PRECISION:: E1, E2, phi1, phi2, E3

	E3=(E1 * phi2 - E2 * phi1) / (phi2 - phi1)
END SUBROUTINE secant

SUBROUTINE tir(nequ, E1, E2, error, a,b,N, yyin)
	IMPLICIT NONE
	DOUBLE PRECISION:: E1, E2, E3, a, b
	DOUBLE PRECISION:: phi1, phi2, phi3, dx, x, error, yyin(nequ), yyout(nequ)
	INTEGER:: i, nequ, N, iter

	i=0 !incialitzo
	iter=500
	!Pas 1 i 2

	do i = 1, iter
		call integrate(nequ, E1, a, b, N,yyin, phi1)
		call integrate(nequ, E2, a, b, N, yyin, phi2)
		print*, phi1, phi2

		call secant(E1, E2, phi1, phi2, E3)
		call integrate(nequ, E3, a, b, N,yyin, phi3)

		if (abs(phi3).LT.error) then
			write(*,*) 'Ha convergit', i, E3
			exit

	    else
	    	E1 = E2
		    E2 = E3
		end if
	end do	

END SUBROUTINE tir

```

Here is my original code, it does compile. Again I have the same error I mentioned.

[Next page](https://fortran-lang.discourse.group/t/help-with-runge-kutta-4/6910.md?page=2)