Adding solutions

.gitignore

@@ -193,3 +193,5 @@ sympy-plots-for-*.tex/
 notebooks/*.ipynb_checkpoints
 notebooks/*dataxy
 notebooks/*gnuxy
+src/*dataxy
+src/*gnuxy

src/solutions/04_newton.f90

+program newton
+
+  implicit none
+
+! Use a Newton iteration to solve the equation
+!
+!  x     -    current approximation to the solution
+!  f(x)  -    polynomial function
+!  df    -    derivative of f with respect to x
+!  xo    -    previous guess for solution
+!  eps   -    convergence criterion
+!  dx    -    change in solution approximation
+!  it    -    number of iterations
+!  itmax -    maximum number of iterations
+
+  integer, parameter :: itmax = 1000 
+  real(8), parameter :: eps   = 1.e-6
+
+  integer :: it
+  real(8) :: x,f,df,xo
+
+  it = 0
+  f  = 1d0
+
+  write(*,*) 'Try to solve "x**3+x-10=0"'
+  write(*,*) 'What is your initial guess for the solution?'
+  read(*,*) x
+
+  do while( dabs(f)>eps .and. it<=itmax )
+    it=it+1
+    xo=x
+    call derivate(xo,f,df)
+    x = xo -f/df
+    write(*,*) 'it = ',it,', x = ',x,', f(x) = ',f
+  end do
+
+end program
+
+! ******************************************************************************************
+
+subroutine derivate(x,f,d)
+
+! Evaluate the function f(x)=x**3+x-10
+! also return the derivative of the function
+                                             
+   implicit none
+
+   real(8), intent(in)  :: x
+   real(8), intent(out) :: f, d
+
+   d = 3*x**2 + 1d0
+   f = x**3 + x - 10d0
+
+end subroutine

src/solutions/04_plot.f90

+program plot
+
+  implicit none
+      
+  real(4) :: fx, x
+  integer :: i
+
+  open (112,file='04_gnuxy')
+
+  write(112,*) 'set grid'
+  write(112,*) 'set xzeroaxis'
+  write(112,*) 'set yzeroaxis'
+  write(112,*) 'set border 0          # remove frame'
+  write(112,*) 'set xtics axis        # place tics on axis rather than on border'
+  write(112,*) 'set ytics axis'
+  write(112,*) 'set ticscale 0        # [optional] labels only, no tics'
+  write(112,*) 'set xtics add ("" 0)  # suppress origin label that lies on top of axis'
+  write(112,*) 'set ytics add ("" 0)  # suppress origin label that lies on top of axis'
+  write(112,*) ''
+  write(112,*) '# if arrows are wanted only in the positive direction'
+  write(112,*) 'set arrow 1 from 0,0 to graph 1, first 0 filled head'
+  write(112,*) 'set arrow 2 from 0,0 to first 0, graph 1 filled head'
+  write(112,*) ''
+  write(112,*) '# if arrows in both directions from the origin are wanted'
+  write(112,*) 'set arrow 3 from 0,0 to graph 0, first 0 filled head'
+  write(112,*) 'set arrow 4 from 0,0 to first 0, graph 0 filled head'
+  write(112,*) ''
+  write(112,*) 'set nokey'
+  write(112,*) 'set xrange [-4:4]'
+  write(112,*) 'plot "04_dataxy_1" using 1:2 with lines lt rgb "blue"'
+  write(112,*) 'pause -1'
+
+  close(112)
+
+  ! Generate x-y pairs for the graph
+
+  open (112,file='04_dataxy_1')
+  do i=-40,40
+    x  = .1*i
+    fx = x**3+x-10.0
+    write(112,*) x, fx
+  end do
+  close(112)
+
+  print *, ' Hit the Return (Enter) key to continue'
+
+  call system ('gnuplot 04_gnuxy')
+  stop
+
+end

src/solutions/09_ChristmasTree.f90

+program ChristmasTree
+
+  implicit none
+
+  integer i, h, hmax, line, ball
+  character(8) hmax_string
+  character(1) sball
+
+  ! get the command line argument
+  call getarg(1,hmax_string)
+
+  ! cast string to integer
+  read(hmax_string,*) hmax
+
+  ball=1 
+  do h=1,hmax
+    line=1
+    ! write spaces to align the head of the tree
+    write(*,'(a)',advance='no') repeat(' ',hmax-h)
+    ! loop to decide when we have to create a new line
+    do while(line.le.(2*h-1))
+      ! modulo to decide when we have to put a ball or not
+      sball='#'
+      if(mod(ball,6).eq.0) sball='o'
+      write(*,'(A)',advance='no') sball
+      ! increment ball/line decision variables
+      ball=ball+1
+      line=line+1
+    enddo
+    ! create a new line
+    write(*,*)
+  enddo
+
+end program

src/solutions/10_eqdiff.f90

+program eqdiff
+
+  implicit none
+
+  integer, parameter :: nt  = 200
+  integer, parameter :: st = 1
+
+  integer :: n
+  real(8) :: pi, A, delta
+
+  real(8), dimension(nt+st) :: U, U1, U2, UU
+  real(8), dimension(nt+st) :: E, E1, E2
+
+  open(100,file="10_results.txt")
+
+  pi=2d0*dacos(0)
+  A=1d0
+  delta=pi/50d0
+
+  U(:)=0d0
+  E(:)=0d0
+  U1(:)=0d0
+  E1(:)=0d0
+  U2(:)=0d0
+  E2(:)=0d0
+  UU(:)=0d0
+
+  write(*,*) "Starting the loop..."
+  write(100,6000) st,U(1),U1(1),U2(1),UU(1),E(1),E1(1),E2(1)
+
+  do n=st,nt
+    U(n+1)=U(n)+delta*A*dcos(n*delta)
+    U1(n+1)=U1(n)+delta*A*dcos((n+0.5d0)*delta)
+    U2(n+1)=U2(n)+delta*A*dcos((n+1d0)*delta)
+
+    UU(n+1)=A*dsin(n*delta)
+
+    E(n+1)=dabs(UU(n+1)-U(n+1))
+    E1(n+1)=dabs(UU(n+1)-U1(n+1))
+    E2(n+1)=dabs(UU(n+1)-U2(n+1))
+
+    write(100,6000) (n+1),U(n+1),U1(n+1),U2(n+1),UU(n+1),E(n+1),E1(n+1),E2(n+1)
+  end do
+
+  write(*,5000) "End of the loop: (",nt,") iterations!"
+
+  open (112,file='10_gnuxy')
+
+  write(112,*) "set multiplot layout 2,1 rowsfirst"
+  write(112,*) "set xlabel 'N'"
+  write(112,*) "set xrange [1:201]"
+  write(112,*)
+  write(112,*) "set ylabel 'U'"
+  write(112,*) "set yrange [-1.2:1.2]"
+  write(112,*) "plot '10_results.txt' using 2 title 'U with N' with lines lt rgb 'blue', \"
+  write(112,*) "     '10_results.txt' using 3 title 'U with N+0.5' with lines lt rgb 'red', \"
+  write(112,*) "     '10_results.txt' using 4 title 'U with N+1' with lines lt rgb 'green', \"
+  write(112,*) "     '10_results.txt' using 5 title 'U analytic' with lines lt rgb 'black'"
+  write(112,*) 
+  write(112,*) "set ylabel 'E'"
+  write(112,*) "set yrange [-0.05:0.2]"
+  write(112,*) "plot '10_results.txt' using 6 title 'E with N' with lines lt rgb 'blue', \"
+  write(112,*) "     '10_results.txt' using 7 title 'E with N+0.5' with lines lt rgb 'red', \"
+  write(112,*) "     '10_results.txt' using 8 title 'E with N+1' with lines lt rgb 'green'"
+  write(112,*) 
+  write(112,*) "pause -1"
+  write(112,*) "unset multiplot"
+
+  write(*,*) "See 10_results.txt"
+  write(*,*)
+  write(*,*) "Hit the Return (Enter) key to continue"
+
+  close(112)
+  close(100)
+
+  call system('gnuplot 10_gnuxy')
+
+  5000 format(1X,A,I3,A)
+  6000 format(1X,I5,7F14.9)
+
+end