title: Matlab differential equation solver ode45 error Index exceeds the number of array elements (1)
body: <p>I'm struggling with solving second order differential equation Matlab. It gives me variety of errors, one of these error logs says</p>
  <p>Index exceeds the number of array elements (1).</p>
  <p>Error in
  symengine&gt;@(x,Y)[Y(2);Y(1)-(Y(2).<em>(pi.</em>(tan(x.<em>(5.0./2.9e+1)).^2.</em>(5.0./2.9e+1)+5.0./2.9e+1).*1.0./sqrt(1.0./Y(1).^2.*8.0+1.0).*5.453697799462048e-1-(1.0./Y(1).^5.*Y(2).*1.821511711048492e-2)./(1.0./Y(1).^2.*8.0+1.0)+1.0./Y(1).^7.*Y(2).*1.0./(1.0./Y(1).^2.*8.0+1.0).^2.*7.286046844193969e-2+pi.<em>tan(x.</em>(5.0./2.9e+1)).*1.0./Y(1).^3.*Y(2).*1.0./(1.0./Y(1).^2.*8.0+1.0).^(3.0./2.0).*4.362958239569638))./((1.0./Y(1).^4.*9.107558555242461e-3)./(1.0./Y(1).^2.*8.0+1.0)+pi.<em>tan(x.</em>(5.0./2.9e+1)).*1.0./sqrt(1.0./Y(1).^2.*8.0+1.0).<em>1.09073955989241-2.0)+(pi.<em>tan(x.</em>(5.0./2.9e+1)).</em>((1.0./Y(1).^2.*1.0./sqrt(1.0./Y(1).^2.<em>8.0+1.0).</em>(1.0./Y(1).^2.*1.438008285978908e+22+7.190041429894541e+21).*5.806568032841007e-25)./(pi.<em>tan(x.</em>(5.0./2.9e+1)))+1.0).*Y(1).^5.*1.0./(1.0./Y(1).^2.*8.0+1.0).^(3.0./2.0).*7.664768928982835e-3)./((1.0./Y(1).^4.*4.55377927762123e-9)./(1.0./Y(1).^2.*8.0+1.0)+pi.<em>tan(x.</em>(5.0./2.9e+1)).*1.0./sqrt(1.0./Y(1).^2.*8.0+1.0).*5.453697799462048e-7-1.0e-6)]</p>
  <p>Error in odearguments (line 90)
  f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.</p>
  <p>Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);</p>
  <p>Error in ode (line 77)
  ySol = ode45(M,interval, 0,1);</p>
  <p>What do I do wrong?</p>
  <pre><code>close all;
  clear;
  
  %constant init, in cgs units
  
  c = 2.99792458*1e10;                                    %speed of light
  %e_0 =1; %8.8541878128*1e-12; SI                       %vacuum permitivitty
  e = 4.8032*1e-10;   % 1.60217663*1e-19 SI;                           % elementary charge
  m = 9.1094*1e-28 ;  %9.10956*1e-31; SI                              %electron rest mass
  
  
  %input init
  part=10000;
  t = linspace(0,1*pi,part);                          % discrete time interval for plotting
  lambda_0 = 1060e-7;     %1060*1e-9;                     %laser wavelength
  r_0 = 10e-4;            %10e-6;                                           % focal spot radius
  r = 0;                                     %radial coordinate
  n_0 = 1.7*1e21;                                            %density in cm^-3
  a_0 = 4;
  % a = a0*exp(-2*t*t/(pi^2));                  % zmen vzorec pro a
  % a = e*abs(E)/(m*w*c);
  
  
  %calculation
  w_0 = 2*pi*c/lambda_0;                                            %frequency of the laser
  % k_0 = (w_0/c)*sqrt(e_0);                                            %wavenumber
  R = w_0*r_0^2/c;
  z = linspace(0,3*R,part);                          % discrete distance interval
  % x = z/R;                                                                    %normalised distance
  % n(1:part) = n_0*tan(x/5.8);                                               %density ramp
  % w_p = sqrt(4*pi*n*e^2/(m));                             %plasma frequency
  
  
  % plot(x,w_p.^2/w_0.^2)
  
  %symbolic calculation
  
  syms x;
  syms Y(x);
  n = n_0*tan(sym(x/5.8));
  w_p = sqrt(4*pi*n*e^2/(m));
  a = sqrt(a_0^2/(Y*Y)*exp(-r^2/(r_0*Y*Y)));
  gamma = sqrt(1+a^2/2);
  e_0 = 1-(w_p^2/w_0^2*(1/(gamma))+c^2*a^2/(w_0^2*r_0^2*Y*Y*gamma^2));
  o = w_p^2/(4*w_0^2)*Y*Y*Y*Y*Y*Y*a^2/gamma^3*(1+(c^2*(8+a^2))/(r_0^2*w_p^2*Y^2*gamma));
  
  eqn = diff(Y,x,2) == 1/Y*Y*Y-1/(2*e_0)*(diff(e_0,x)*diff(Y,x))-R^2/(r_0^2*e_0)*o*Y;
  
  dY=diff(Y);
  cond1 = Y(0) == 1;
  cond2 = dY(0) == 0;
  conds = [cond1 cond2];
  
  interval = [0 3*R];
  
  
  [V] = odeToVectorField(eqn);
   
  M = matlabFunction(V,'vars',{'x','Y'});
  t0=0;
  y0=1;
  yp0=0;
  ySol = ode45(M,interval, 0,1);
  </code></pre>