title: Algorithm 11.2 Nonlinear Shooting Method (Burden and Faires) Python
body: <p>I am trying to program nonlinear shooting method based on algorithm 11.2 from Numerical Analysis (Burden and Faires). However, after running the program, the numerical result i am getting is different from the answer in the textbook. i think there is something wrong in my coding but i cannot figure it out. I attached actual algorithm in the picture. <a href="https://i.stack.imgur.com/zfL2K.png" rel="nofollow noreferrer">Algorithm 11.2</a>
  <a href="https://i.stack.imgur.com/4AbPt.png" rel="nofollow noreferrer">Algorithm 11.2</a> 
  <a href="https://i.stack.imgur.com/7Ypbg.png" rel="nofollow noreferrer">Algorithm 11.2</a></p>
  
  <p>Here is the code</p>
  
  <pre><code>from numpy import zeros, abs
  
  def shoot_nonlinear(a,b,alpha, beta, n, tol, M):
  
      w1 = zeros(n+1)  
      w2 = zeros(n+1)
      h = (b-a)/n
      k = 1
      TK = (beta - alpha)/(b - a)
  
      print("i""  x" "     " "W1""     " "W2")
      while k &lt;= M:
  
          w1[0] = alpha
          w2[0] = TK
          u1    = 0
          u2    = 1
  
          for i in range(1,n+1):
              x = a + (i-1)*h     #step 5
  
              t = x + 0.5*(h)
  
              k11 = h*w2[i-1]     #step 6
  
              k12 = h*f(x,w1[i-1],w2[i-1])
              k21 = h*(w2[i-1] + (1/2)*k12)
              k22 = h*f(t, w1[i-1] + (1/2)*k11, w2[i-1] + (1/2)*k12)
              k31 = h*(w2[i-1] + (1/2)*k22)
              k32 = h*f(t, w1[i-1] + (1/2)*k21, w2[i-1] + (1/2)*k22)
              t   = x + h
              k41 = h*(w2[i-1]+k32)
              k42 = h*f(t, w1[i-1] + k31, w2[i-1] + k32)
              w1[i] = w1[i-1] + (k11 + 2*k21 + 2*k31 + k41)/6
              w2[i] = w2[i-1] + (k12 + 2*k22 + 2*k32 + k42)/6   
              kp11 = h*u2
              kp12 = h*(fy(x,w1[i-1],w2[i-1])*u1 + fyp(x,w1[i-1], w2[i-1])*u2)
              t    = x + 0.5*(h)
              kp21 = h*(u2 + (1/2)*kp12)
              kp22 = h*((fy(t, w1[i-1],w2[i-1])*(u1 + (1/2)*kp11)) + fyp(x+h/2, w1[i-1],w2[i-1])*(u2 +(1/2)*kp12))
              kp31 = h*(u2 + (1/2)*kp22)
              kp32 = h*((fy(t, w1[i-1],w2[i-1])*(u1 + (1/2)*kp21)) + fyp(x+h/2, w1[i-1],w2[i-1])*(u2 +(1/2)*kp22))
              t    = x + h
              kp41 = h*(u2 + kp32)
              kp42 = h*(fy(t, w1[i-1], w2[i-1])*(u1+kp31) + fyp(x + h, w1[i-1], w2[i-1])*(u2 + kp32))
              u1 = u1 + (1/6)*(kp11 + 2*kp21 + 2*kp31 + kp41)
              u2 = u2 + (1/6)*(kp12 + 2*kp22 + 2*kp32 + kp42)
  
  
          r = abs(w1[n]) - beta
          #print(r)
          if r &lt; tol:
              for i in range(0,n+1):
                  x = a + i*h
                  print("%.2f %.2f %.4f %.4f" %(i,x,w1[i],w2[i]))
              return
  
          TK = TK -(w1[n]-beta)/u1
  
          k = k+1
  
  
      print("Maximum number of iterations exceeded")   
      return
  </code></pre>
  
  <h1>function for 2nd order boundary value problem</h1>
  
  <pre><code>def f(x,y,yp):
  fx = (1/8)*(32 + 2*x**3 -y*yp)
  return fx
  
  def fy(xp,z,zp):
  fyy = -(1/8)*(zp)
  return fyy
  
  def fyp(xpp,zpp,zppp):
  fypp = -(1/8)*(zpp)
  return fypp
  
  a = 1         # start point
  b = 3         # end point
  alpha = 17    # boundary condition
  beta = 43/3   # boundary condition
  N = 20        # number of subintervals
  M = 10        # maximum number of iterations
  tol = 0.00001 # tolerance
  
  
  shoot_nonlinear(a,b,alpha,beta,N,tol,M)
  </code></pre>
  
  <p>My result</p>
  
  <pre><code>i      x     W1      W2
  0.00 1.00 17.0000 -16.2058
  1.00 1.10 15.5557 -12.8379
  2.00 1.20 14.4067 -10.2482
  3.00 1.30 13.4882 -8.1979
  4.00 1.40 12.7544 -6.5327
  5.00 1.50 12.1723 -5.1496
  6.00 1.60 11.7175 -3.9773
  7.00 1.70 11.3715 -2.9656
  8.00 1.80 11.1203 -2.0783
  9.00 1.90 10.9526 -1.2886
  10.00 2.00 10.8600 -0.5768
  11.00 2.10 10.8352 0.0723
  12.00 2.20 10.8727 0.6700
  13.00 2.30 10.9678 1.2251
  14.00 2.40 11.1165 1.7444
  15.00 2.50 11.3157 2.2331
  16.00 2.60 11.5623 2.6951
  17.00 2.70 11.8539 3.1337
  18.00 2.80 12.1883 3.5513
  19.00 2.90 12.5635 3.9498
  20.00 3.00 12.9777 4.3306
  </code></pre>
  
  <p>Actual result for w1</p>
  
  <pre><code>x     W1
  1.0   17.0000
  1.1   15.7555
  1.2   14.7734
  1.3   13.3886
  1.4   12.9167
  1.5   12.5601
  1.6   12.3018
  1.7   12.1289
  1.8   12.0311
  1.9   12.0000
  2.0   12.0291
  2.1   12.1127
  2.2   12.2465
  2.3   12.4267 
  2.4   12.6500
  2.5   12.9139
  2.6   13.2159
  2.7   13.5543
  2.8   13.9272
  2.9   14.3333
  3.0   14.7713
  </code></pre>