Analytical Jacobian for BVP4C. How to write it ?

This is my code for the non-linear differential equation as given below in the image :-

function [sol] = pbeexp()
  el = 4.8*10^-10;  %electron charge in StatCoulomb in CGS
  eps = 80; %permittivity of water at 20 degrees C
  k = 1.3807*10^-16;   %Boltzmann constant in CGS
  T = 293 ;          %absolute temperature at 20 degrees C
  l = el^2/(eps*k*T);%Bjerrum length in cm
  pn = k*T/el; %non dimensional potential
  psiwall =300*10^-3*0.00333564; %Electric field at the electrode in StatVolts
  %Formulating the problem
  x=linspace(0,35*10^-7,10000)/l;%Linearly distributed grid
  solinit = bvpinit(x,[psiwall/pn 0]); % Initial mesh
  options = bvpset('RelTol',1e-5,'AbsTol',1e-5); % Tolerance values
  sol = bvp4c(@pbeode,@pbebc,solinit,options);%Initial solution structure
  xint=linspace(sol.x(1),sol.x(end),size(sol.x,2));%Linearly distributed grid
  sxint=deval(sol,xint);
  dist = xint;
  pot = sxint(1,:);
  plot(dist,pot)
  function dydx = pbeode(x,psi)
      el = 4.8*10^-10;  %electron charge in StatCoulomb in CGS
      eps = 80; %permittivity of water at 20 degrees C
      k = 1.3807*10^-16;   %Boltzmann constant in CGS
      T = 293 ;          %absolute temperature at 20 degrees C
      l = el^2/(eps*k*T);%Bjerrum length in cm
      H = 7.8*10^-7;      %length of probe in cm
      a3 = (1/4.6)*10^3 ;          %steric packing in cc
      czero = 6.023*10^20*0.1 ;       %bulk counterion concentration
      nu = 2*czero*a3;    %steric size parameter for use in solution
      %% Non - dimensional quantities
      pn = k*T/el; %non dimensional potential
      %% PBE terms for different components of the model
      A1 = (4*pi*czero*l*sinh(psi(1)))/(1+2*nu*((sinh(psi(1)/2)^2)));% expression for counterions in the solution 
      dydx = [ psi(2)
          (2*l^2*A1)];
  end 
  %% Providing the B.C.s
  function bc = pbebc(psia,psib)
      psiwall = 300*10^-3*0.00333564; %Electric field at the electrode in StatVolts
      pn = k*T/el; %non dimensional potential
      bc = [ psia(1)-psiwall/pn
          psib(1)];
  end
  end

*I wish to provide "analytical Jacobian" for BVP4C solver. Can anybody tell me how to do it ? Any hints or suggestions is available. I checked the documentation but I could not write the Jacobian for my Non Linear ODE.

Thanks.*