### A sample of toy FDM code

Solving equation

$u_t-u''=f(x,u)$

$f(x,u)=-2+2*(\sin(u)-\sin(x^2))$

Time integration is full implicit, and during each time step Newton’s method is employed.

PROGRAM MAIN
!--------------------------------------------------
!Implicit time evolution of equation
!u_t-u''=f(x,u)
!f(x,u)=-2+r*(sin(u)-sin(x^2)),r=2
!--------------------------------------------------
! t=time points
! dt=time step
! h=neighboring nodes distance
! u=unknown
! x=node location
! n=No. of unknowns

!USE SUNPERF

IMPLICIT NONE
INTEGER(4) :: n
REAL(4) :: elapse,time(2),etime
REAL(8) :: t=0,dt=0.1
REAL(8) :: h
REAL(8),ALLOCATABLE :: u(:),x(:)
INTEGER(4) :: count,i

n=50
h=1/(DBLE(n)+1.)
ALLOCATE (x(n),u(n))
x=(/(i*h,i=1,n)/)
u=(/((i*h)**0.5,i=1,n)/)

OPEN(UNIT=100,FILE='data',STATUS='UNKNOWN',ACTION='&
&WRITE')
DO i=1,n
WRITE(100,*) t,x(i),u(i)
END DO

DO WHILE (t<1.0)
CALL newton(n,h,dt,u,x)
t=t+dt
DO i=1,n
WRITE(100,*) t,x(i),u(i)
END DO END DO CLOSE(UNIT=100)
elapse=etime(time)
WRITE(*,*) 'elapsed:',elapse
END PROGRAM MAIN
!--------------------------------------------------
SUBROUTINE newton(n,h,dt,u,x)
IMPLICIT NONE
INTEGER(4),INTENT(IN) :: n
REAL(8),INTENT(IN) :: x(n),h,dt
REAL(8),INTENT(INOUT) :: u(n)
REAL(8) :: r=2,c(n),error,b(n,1),b0(n,1),a(2,n)
REAL(8) :: alpha=1.0,beta=-1.0
INTEGER(4) :: info,i,iter
c=u*h*h/dt
iter=0
error=1
DO WHILE (error>1.e-13)
iter=iter+1
a(1,:)=-1
a(2,:)=2+h*h/dt
b(:,1)=(-2+r*(SIN(u)-SIN(x*x)))*h*h+c
b(n,1)=b(n,1)+1
CALL DSBMV('u',n,1,alpha,a,2,u,1,beta,b,1)
b=-b
a(2,:)=a(2,:)-r*COS(u)*h*h
CALL DGTSV(n,1,a(1,2:),a(2,:),a(1,2:),b,n,info)
error=MAXVAL(ABS(b))
u=u+b(:,1)
END DO
WRITE(*,*) 'INFO',info,'Number of iterations:',iter
END SUBROUTINE newton
!--------------------------------------------------


Output would be like:

% f95 -dalign implicit_newton.f95 -xlic_lib=sunperf
% a.out
INFO 0 Number of iterations: 4
INFO 0 Number of iterations: 4
INFO 0 Number of iterations: 4
INFO 0 Number of iterations: 4
INFO 0 Number of iterations: 4
INFO 0 Number of iterations: 3
INFO 0 Number of iterations: 3
INFO 0 Number of iterations: 3
INFO 0 Number of iterations: 3
INFO 0 Number of iterations: 3
elapsed: 0.061407


Of course the exact solution is $u(x)=x^2$. Here is time evolution plot of the numerical solution.