Create runge_kutta module, begin with simple 2-stage routines

Makefile

@@ -18,11 +18,14 @@ FLIBS = -lblas -llapack
 .DEFAULT_GOAL := $(BIN)/main
 
 # Dependencies of main program
-objects=$(OBJ)/misc.o
+objects=$(OBJ)/misc.o \
+	$(OBJ)/runge_kutta.o
 
 # Modules
 $(OBJ)/misc.o: $(SRC)/misc.f90
 	$(FF) $(FFLAGS) -J$(OBJ) -c -o $@ $<
+$(OBJ)/runge_kutta.o: $(SRC)/runge_kutta.f90
+	$(FF) $(FFLAGS) -J$(OBJ) -c -o $@ $<
 
 # Main program
 $(OBJ)/main.o: $(SRC)/main.f90 $(objects)

src/main.f90

@@ -1,16 +1,28 @@
 program main
   use iso_fortran_env, only: wp => real64
   use misc, only: myfun
+  use runge_kutta, only: midpoint, heun, ralston
   implicit none
 
   integer :: ii
-  real(wp) :: t, y, dy
 
-  t = 0_wp
-  y = 0_wp
+  ! real(wp) :: t, y, dy
+  ! t = 0_wp
+  ! y = 0_wp
+  ! dy = myfun(t, y)
+  ! print *, t, y, dy
 
-  dy = myfun(t, y)
+  integer, parameter :: N = 2
+  real(wp), dimension(N) :: c, b
+  real(wp), dimension(N,N) :: a
 
-  print *, t, y, dy
+  call ralston(N, a, b, c)
+
+  print *, c
+  print *, b
+
+  do ii = 1, N
+    print *, a(ii, :)
+  end do
 
 end program main

src/runge_kutta.f90

@@ -0,0 +1,54 @@
+module runge_kutta
+  use iso_fortran_env, only: wp => real64
+  implicit none
+
+contains
+
+  subroutine midpoint(N, a, b, c)
+    integer, intent(in) :: N
+    real(wp) :: alpha
+    real(wp), intent(out), dimension(N) :: c, b
+    real(wp), intent(out), dimension(N,N) :: a
+
+    alpha = 0.5_wp
+    call two_stage(N, a, b, c, alpha)
+
+  end subroutine midpoint
+
+  subroutine heun(N, a, b, c)
+    integer, intent(in) :: N
+    real(wp) :: alpha
+    real(wp), intent(out), dimension(N) :: c, b
+    real(wp), intent(out), dimension(N,N) :: a
+
+    alpha = 1._wp
+    call two_stage(N, a, b, c, alpha)
+
+  end subroutine heun
+
+  subroutine ralston(N, a, b, c)
+    integer, intent(in) :: N
+    real(wp) :: alpha
+    real(wp), intent(out), dimension(N) :: c, b
+    real(wp), intent(out), dimension(N,N) :: a
+
+    alpha = 2._wp/3._wp
+    call two_stage(N, a, b, c, alpha)
+
+  end subroutine ralston
+
+  subroutine two_stage(N, a, b, c, alpha)
+    integer, intent(in) :: N
+    real(wp), intent(in) :: alpha
+    real(wp), intent(out), dimension(N) :: c, b
+    real(wp), intent(out), dimension(N,N) :: a
+
+    c = [0._wp, alpha]
+    b = [(1._wp - (1._wp/(2._wp*alpha))), (1._wp/(2._wp*alpha))]
+
+    a = 0._wp
+    a(2,1) = alpha
+
+  end subroutine two_stage
+
+end module runge_kutta