Plotting the results of the differential equation plotted suing RK4 method

Question

I have a code for the solution of the Lane-Emden equation which is giving correct values after solving but when I plot those lists using matplotlib, I am getting blank images and not correct plots.

import numpy as np
import matplotlib.pyplot as plt

n = 14
theta_0 = 1
phi_0 = 0
h = 0.01
xi_0 = 0
xi_max = 100
   
theta = theta_0
phi = phi_0
xi = xi_0 + h
   
Theta = [[] for i in range(n)]
Phi = [[] for i in range(n)]
Xi = [[] for i in range(n)]
    
for i in range(n):
    Theta[i].append(theta)
    Phi[i].append(phi)
    Xi[i].append(xi)
     
def dThetadXi(phi,xi): #r1
    return -phi/xi**2
        
def r2(phi,xi):
    return dThetadXi(phi+h,xi+h*dThetadXi(phi,xi))
    
def r3(phi,xi):
    return dThetadXi(phi+h,xi+h*r2(phi,xi))
      
def r4(phi,xi):
    return dThetadXi(phi+h,xi+h*r3(phi,xi))
    
def dPhidXi(theta,xi,n): #k1
    return theta**(n)*(xi**2)
    
        
def k2(theta,xi,n):
    return dPhidXi(theta+h,xi+h*dPhidXi(theta,xi,n),n)
    
def k3(theta,xi,n):
    return dPhidXi(theta+h,xi+h*k2(theta,xi,n),n)
        
def k4(theta,xi,n):
    return dPhidXi(theta+h,xi+h*k3(theta,xi,n),n)
        
for i in range(n):  
    while xi < xi_max:
        if theta < 0:
            break
        dTheta = (step/6)*(dThetadXi(phi,xi)+2*r2(phi,xi)+2*r3(phi,xi)+r4(phi,xi))
        dPhi = (step/6)*(dPhidXi(theta,xi,i/2.)+2*k2(theta,xi,n)+2*k3(theta,xi,n)+k4(theta,xi,n))
        theta = theta+ dTheta
        phi = phi +dPhi
        xi = xi + h
        Theta[i].append(theta)
        Phi[i].append(phi)
        Xi[i].append(xi)
    print (i/2., round(xi,2), round(dThetadXi(phi,xi),2), round(xi/3./dThetadXi(phi,xi),2), round(1./(4*np.pi*(i/2.+1))/dThetadXi(phi,xi)**2,2))
    theta = theta_0
    phi = phi_0
    xi = xi_0 + h
    plt.plot(phi, xi)
    plt.show()

The correct plots should like this. Thanks.

EDIT After implementing the edits suggested below, I am getting this error output.

usr/local/lib/python3.7/dist-packages/numpy/core/_asarray.py:136: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  return array(a, dtype, copy=False, order=order, subok=True)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
TypeError: float() argument must be a string or a number, not 'list'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
<ipython-input-20-df5bae47b875> in <module>()
     62     phi = phi0
     63     xi = xi0 + step
---> 64     plt.plot(Phi, Xi)
     65 

7 frames
/usr/local/lib/python3.7/dist-packages/numpy/core/_asarray.py in asarray(a, dtype, order)
     81 
     82     """
---> 83     return array(a, dtype, copy=False, order=order)
     84 
     85 

ValueError: setting an array element with a sequence.

EDIT2

Answer 1

To summarize the comments under the original post which led to the solution:

Consider calling plt.plot() and plt.show() only once after the for-loop and plot individual sublists of Xi and Phi with plt.plot(Xi[0],Phi[0]) instead of single points like you did here plt.plot(xi, phi).