title: Python: use odeint to solve second-order linear differential equation
body: <p>I have two differential equations like this: </p>
  
  <ol>
  <li><code>al=arctan(ri')</code></li>
  <li><code>ri''=(vz-C_alz(ri,al))*(1+ri' **2)/C_aln(ri,al)</code></li>
  </ol>
  
  <p>The initial value is: <code>ri(l=0)=2.5, ri'(l=0)=0</code></p>
  
  <p><code>vz</code> is a function about l and <code>C_alz</code>, <code>C_aln</code> are two functions about <code>ri</code> and <code>al</code>. These 3 complex mathematical expressions are calculated through the code 1, the results are correct and the trajectory of <code>ri</code> can be used as the reference to the result of code 2 .</p>
  
  <pre><code># -*- coding: utf-8 -*-
  """
  Created on Fri Dec 18 14:41:07 2015
  @author: marswang
  code 1
  """
  
  from numpy.linalg import inv
  from sympy import sqrt, sin, cos, tan, atan, diff
  from numpy import array
  import sympy as sp
  import matplotlib.pyplot as pl
  import numpy as np
  def r_trajektor (L,R):  
      B= np.array([[1, L[0], L[0]**2,L[0]**3 ], [0, 1, 2*L[0], 3*L[0]**2], [1, L[1], L[1]**2, L[1]**3], [0, 1, 2*L[1], 3*L[1]**2]])
      X= np.array([[R[0]], [0], [R[1]], [0]])
      A=np.dot(inv(B), X)
      return A    
  
  a=7.6
  rh_s=2420                          
  rh_m=2285                          
  the_0=11                           
  Rc=49                              
  vc= 0                             
  vp=10
  
  # two points in the r-trajektory and 1. derivative:
  # point1 l0=0,r(l0)=2.5 and r'(l0)=0 
  # point2 l1=36,r(l1)=22.5 and r'(l1)=0
  L= array([[0], [36]])
  R= array([[2.5], [22.5]])
  # A symbolic derivation expression of vg
  ri=sp.Symbol('ri') 
  al=sp.Symbol('al')
  l=sp.Symbol('l')  
  A=r_trajektor(L,R)
  # convert 'list' to 'float'
  rir=float(A[0])+float(A[1])*l+float(A[2])*l**2+float(A[3])*l**3
  dri=diff(rir,l)
  alr=atan(diff(rir,l))
  dal=diff(alr,l)
  hi=a*sqrt((1-sin(the_0+al))/(1+a/(sqrt(2)*ri)))
  hir=hi.subs([(al,alr),(ri,rir)])
  C_aln=(1-ri**2/(Rc**2))*diff(hi,al)+(a**2)*ri/(Rc**2)*sin(the_0+al)
  print C_aln
  C_alz=1-(rh_s*ri**2)/(rh_m*Rc**2)+((1-(ri**2)/(Rc**2))*diff(hi,ri)-2*ri*hi/(Rc**2)-(a**2)/(Rc**2)*cos(the_0+al))*tan(al)
  print C_alz
  vg=((vp-vc)/(diff(al,l)*C_aln+C_alz)).subs([(al,alr),(ri,rir)])
  vz=10/vg
  print vz
  
  
  # the function has to be worked with the arrays from numpy too.
  # Thus have to lambdify the expression
  
  l_val = np.linspace(0.0,36.0,100.0)
  
  fr=sp.lambdify(l,rir,"numpy")
  r_val=fr(l_val)
  ## plot
  pl.figure(figsize=(10,10),dpi=98)
  p1 =pl.subplot(2, 1, 1)
  p1.plot(l_val,r_val,"g-",label="ri",)
  p1.axis([0.0, 40.0,0.0,25.0])
  p1.set_ylabel("ri in mm",fontsize=14)
  p1.set_xlabel("l in mm",fontsize=14)
  p1.legend(loc=4)
  </code></pre>
  
  <p>Now I want to solve these two differential equations and use code 2 to get the the graphic (trajectory) of <code>ri</code> in relation to <code>l</code> from the calculated <code>vz</code> of code 1 reversely.</p>
  
  <pre><code># -*- coding: utf-8 -*-
  """
  Created on Sat Jan 02 13:15:19 2016
  @author: marswang
  code 2
  """
  from sympy import sqrt, sin, cos, tan, atan
  import matplotlib.pyplot as pl
  import numpy as np
  from scipy.integrate import odeint
  a=7.6                              
  rh_s=2420                          
  rh_m=2285                          
  the_0=11                           
  Rc=49                              
  vc= 0                             
  vp=10
  
  def func(ri,l):
      vz=((-0.00257201646090535*l**2 + 0.0925925925925926*l + 9.25185853854297e-17)*(-0.00633069554352353*sqrt((-sin(atan(-0.00257201646090535*l**2 + 0.0925925925925926*l + 9.25185853854297e-17) + 11) + 1)/(1 + 3.8*sqrt(2)/(-0.000857338820301783*l**3 + 0.0462962962962963*l**2 + 9.25185853854297e-17*l + 2.5)))*(-0.000857338820301783*l**3 + 0.0462962962962963*l**2 + 9.25185853854297e-17*l + 2.5) + 14.44*sqrt(2)*sqrt((-sin(atan(-0.00257201646090535*l**2 + 0.0925925925925926*l + 9.25185853854297e-17) + 11) + 1)/(1 + 3.8*sqrt(2)/(-0.000857338820301783*l**3 + 0.0462962962962963*l**2 + 9.25185853854297e-17*l + 2.5)))*(-(-0.000857338820301783*l**3 + 0.0462962962962963*l**2 + 9.25185853854297e-17*l + 2.5)**2/2401 + 1)/((1 + 3.8*sqrt(2)/(-0.000857338820301783*l**3 +
      0.0462962962962963*l**2 + 9.25185853854297e-17*l + 2.5))*(-0.000857338820301783*l**3 + 0.0462962962962963*l**2 + 9.25185853854297e-17*l + 2.5)**2) - 0.0240566430653894*cos(atan(-0.00257201646090535*l**2 + 0.0925925925925926*l + 9.25185853854297e-17) + 11)) - 484*(-0.000857338820301783*l**3 + 0.0462962962962963*l**2 +
      9.25185853854297e-17*l + 2.5)**2/1097257 + 1)    
      ri0=ri[0]
      ri1=ri[1]       
      C_aln=0.0240566430653894*ri0*sin(atan(ri1) + 11) - 3.8*sqrt((-sin(atan(ri1) + 11) + 1)/(1 + 3.8*sqrt(2)/ri0))*(-ri0**2/2401 + 1)*cos(atan(ri1) + 11)/(-sin(atan(ri1) + 11) + 1)
      C_alz=-484*ri0**2/1097257 + (-0.00633069554352353*ri0*sqrt((-sin(atan(ri1) + 11) + 1)/(1 +
      3.8*sqrt(2)/ri0)) - 0.0240566430653894*cos(atan(ri1) + 11) + 14.44*sqrt(2)*sqrt((-sin(atan(ri1) + 11) + 1)/(1 + 3.8*sqrt(2)/ri0))*(-ri0**2/2401 + 1)/(ri0**2*(1 + 3.8*sqrt(2)/ri0)))*tan(atan(ri1)) + 1    
      ri2=(vz-C_alz)*(1+ri1**2)/C_aln/
      return [ri1, ri2]
  
  init = np.array([2.5, 0.1])
  l = np.linspace(0.0,36.0,100.0)
  sol=odeint(func, init, l)
  pl.plot(l, sol[:,0], color='b')
  pl.legend()
  pl.xlabel(" l in mm ", fontsize=14)
  pl.ylabel(" ri in mm ", fontsize=14)
  pl.show()
  </code></pre>
  
  <p>I have checked the expression of <code>vz</code>, <code>C_aln</code>, <code>C_alz</code> several times. Those expressions come directly from code 1. But result of code 2 is incorrect and what I get is always different from the reference. 
  The reference result of <code>ri</code> should be like the first graphic:</p>
  
  <p><a href="https://i.stack.imgur.com/NplNv.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/NplNv.png" alt="enter image description here"></a></p>
  
  <p>but the calculated result of 'ri' from code 2 is like this:</p>
  
  <p><a href="https://i.stack.imgur.com/HCBLs.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/HCBLs.png" alt="enter image description here"></a></p>
  
  <p>I am not sure if I had used <code>odeint</code> correctly in this situation, what is the reason for this error between these two graphics?</p>
  
  <p>I am new to Python world, any help would be quite appreciated. Thank you.</p>