MATLAB: What is the meaning of “Index in position 1 exceeds array bounds (must not exceed 1)”.

integro-differential equationode 15s

I was trying to solve an integrodifferential equation using MATLAB. I decided to use symbolic to solve ODE, so I wrote the following code. 
%% EQUATIONS GOVERNING THE BUBBLE DYNAMICS
% Yihan Zhang
clc;clear;close all
syms t r R(t)
%% parameters
tau_y=8.6; % Pa



gamma=0.07; % Pa m
rho=1000; % kg/m^3
G=81.5; % Pa
P0=1.13e5; % Pa
deltaP=100; % Pa
eta_s=0.001; % Pa s
R0=100e-4; % m（boundary condition) R(t=0)
Rp0=0; % m/s (boundary condition) dR/dt(t=0)
w=((3*P0+4*gamma/R0+4*G)/(rho*R0^2))^0.5; % s^-1
%% expressions of other variables
tau_rr=-4*G.*(R(t).^3-R0^3)./(3.*r.^3);
tau_tt=-tau_rr./2;
Pg=(P0+2*gamma/R0).*(R0./R(t)).^3;
P_inf=P0+deltaP.*sin(w.*t);
P=Pg+tau_rr-eta_s.*4.*R(t)./R(t)-2.*gamma./R(t);
%% differential equation
eqn= [1000*(diff(R(t),t,2)*R(t)+3/2*diff(R(t),t)^2)==P-P_inf-tau_rr+2*int((tau_rr-tau_tt)/r,R(t),9999)]
[eqs,vars]=reduceDifferentialOrder(eqn,R(t))
[M,F]=massMatrixForm(eqs,vars)
M=odeFunction(M,vars)
Fun=odeFunction(F,vars,r)
H=@(t,R) Fun(t,R)
opt=odeset('mass',M,'InitialSlope',Rp0);
ode15s(H,[0,0.0001],R0,opt)
When I ran this code, the error message is followed:
Index in position 1 exceeds array bounds (must not exceed 1).
Error in
symengine>@(t,in2,param2)[sin(t.*1.842156345156404e3).*-1.0e2-in2(2,:).^2.*1.5e3-(7.0./5.0e1)./in2(1,:)+1.0./in2(1,:).^3.*1.13014e-1+1.0./param2.^4.*(in2(1,:).^3.*3.26e2-3.26e-4).*(in2(1,:)-9.999e3)-1.13000004e5;-in2(2,:)]
Error in Bubble_dynamics2>@(t,R)Fun(t,R)
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
Error in ode15s (line 150)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in Bubble_dynamics2 (line 34)
ode15s(H,[0,0.0001],R0,opt)
I have checked the size of H at line 34, but it is 1*1 sym, which does not exceed 1. I am very confused about the error, but I need to find a way to solve it.
Thank you!

Best Answer

Use the flow shown in https://www.mathworks.com/help/symbolic/odefunction.html#bvlo6zg-1
Replace
[M,F]=massMatrixForm(eqs,vars)
M=odeFunction(M,vars)
Fun=odeFunction(F,vars,r)
H=@(t,R) Fun(t,R)
with
[M,F]=massMatrixForm(eqs,vars)
f = M\F;
H = odeFunction(f,vars)
However this will fail because you have the unresolved variable "r" that is distinct from R(t) . Also your equations involve R(t) and DR(t) so you need at least two boundary conditions but your R0 is only a single boundary condition.
Related Question