I don't know why I am getting the error after all the hard work

% patched-conic gravity assist interplanetary
% trajectory design and optimization

% JPL ephemeris and 64 bit SNOPT algorithm

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;

global emu smu xmu aunit ip1 ip2 ip3 jdtdb0 jdtdbi1

global iephem km ephname eq2000 pmu rep rsoi fbalt

global dv_launch dvm_launch dv_arrival dvm_arrival vinfm_in vinfm_out

global fbrp c3launch jdtdb1 jdtdb2

global otype rito1 vito1 rito2 vito2 dvh

global vinf_in vinf_out jdtdb_soi rp2sc vp2sc

% J2000 ecliptic-to-equatorial transformation matrix

eq2000 = [[1.000000000000000 0  0]; ...
    [0   0.917482062069182  -0.397777155931914]; ...
    [0   0.397777155931914   0.917482062069182]];

% define Astronomical Unit (kilometers)

aunit = 149597870.691;

% radians to degrees conversion factor

rtd = 180.0 / pi;

% define "reference" tdb julian date

jdtdb0 = julian(1, 1, 2000);

% initialize jpl ephemeris

ephname = 'de421.bin';

iephem = 1;

km = 1;

% define planet name vector

pdata = ['Mercury'; 'Venus  '; 'Earth  '; 'Mars   '; ...
    'Jupiter'; 'Saturn '; 'Uranus '; 'Neptune'; 'Pluto  '];

pname = cellstr(pdata);

% planet gravitational constants (kilometer^3/second^2)

pmu(1) = 22032.08;
pmu(2) = 324858.599;
pmu(3) = 398600.436;
pmu(4) = 42828.314;
pmu(5) = 126712767.863;
pmu(6) = 37940626.063;
pmu(7) = 5794549.007;
pmu(8) = 6836534.064;
pmu(9) = 981.601;

xmu = 0.295912208285591149e-03;

smu = xmu * aunit^3 / 86400.0^2;

xmu = 0.899701152970881167e-09;

emu = xmu * aunit^3 / 86400.0^2;

% planet equatorial radius (kilometers)

rep(1) = 2439.7;
rep(2) = 6051.9;
rep(3) = 6378.14;
rep(4) = 3397.0;
rep(5) = 71492.0;
rep(6) = 60268.0;
rep(7) = 25559.0;
rep(8) = 24764.0;
rep(9) = 1151.0;

% sphere-of-influence for each planet (kilometers)

rsoi(1) = 112409.906759936;
rsoi(2) = 616277.129297850;
rsoi(3) = 924647.107795632;
rsoi(4) = 577231.618115568;
rsoi(5) = 48204698.6886979;
rsoi(6) = 54553723.6086354;
rsoi(7) = 51771106.3741412;
rsoi(8) = 86696170.2674129;
rsoi(9) = 15110628.1503479;

% begin simulation

clc; home;

fprintf('\nprogram flyby_matlab\n');

fprintf('\npatched-conic gravity-assist trajectory analysis\n\n');

% request departure tdb calendar date

fprintf('\ndeparture - calendar date guess\n');

('1 <= month <= 12, 1 <= day <= 31, year = all digits!');
3; 5; 2030;
%['\03, 05, 2030\n'] = getdate;

% compute departure tdb julian date

jdtdbi1 = julian(5, 5, 2030);

while(1)
    
    fprintf('\nplease input the departure date search boundary in days\n');
    
    ddays1 = input('30');
    
    if (ddays1 >= 0.0)
        break;
    end
    
end

% request flyby tdb calendar date

fprintf('\n\nflyby - calendar date guess\n');

('1 <= month <= 12, 1 <= day <= 31, year = all digits!');
3; 5; 2037;
%[month, day, year] = getdate;

% compute flyby tdb julian date
jdtdbi2 = julian(03,05,2037);
%jdtdbi2 = julian(month, day, year);

while(1)
    
    fprintf('\nplease input the flyby date search boundary in days\n');
    
    ddays2 = input('?');
    
    if (ddays2 >= 0.0)
        break;
    end
    
end

% request arrival tdb calendar date

fprintf('\n\narrival - calendar date guess\n');

('1 <= month <= 12, 1 <= day <= 31, year = all digits!');
8; 12; 2050;
%(month; day; year) = getdate;

% compute arrival tdb julian date
jdtdbi3 = julian(8, 12, 2050);
%jdtdbi3 = julian(month, day, year);

while(1)
    
    fprintf('\nplease input the arrival date search boundary in days\n');
    
    ddays3 = input('?');
    
    if (ddays3 >= 0.0)
        break;
    end
    
end

% request departure, flyby and arrival planets

for i = 1:1:3
    
    fprintf('\n planet menu\n');
    
    fprintf('\n  <1> Mercury');
    fprintf('\n  <2> Venus');
    fprintf('\n  <3> Earth');
    fprintf('\n  <4> Mars');
    fprintf('\n  <5> Jupiter');
    fprintf('\n  <6> Saturn');
    fprintf('\n  <7> Uranus');
    fprintf('\n  <8> Neptune');
    fprintf('\n  <9> Pluto');
    
    if (i == 1)
        
        while(1)
            
            fprintf('\n\nplease select the departure planet\n');
            
            ip1 = input('?');
            
            if (ip1 >= 1 && ip1 <= 9)
                break;
            end
            
        end
        
    end
    
    if (i == 2)
        
        while(1)
            
            fprintf('\n\nplease select the flyby planet\n');
            
            ip2 = input('?');
            
            if (ip2 >= 1 && ip2 <= 9)
                break;
            end
            
        end
        
    end
    
    if (i == 3)
        
        while(1)
            
            fprintf('\n\nplease select the arrival planet\n');
            
            ip3 = input('?');
            
            if (ip3 >= 1 && ip3 <= 9)
                break;
            end
            
        end
        
    end
    
end

% request lower and upper bounds for the flyby altitude (kilometers)

while(1)
    
    fprintf('\n\nplease input the lower bound for the flyby altitude (kilometers)\n');
    
    fbalt_lwr = input('?');
    
    if (fbalt_lwr >= 0.0)
        
        break;
        
    end
    
end

while(1)
    
    fprintf('\n\nplease input the upper bound for the flyby altitude (kilometers)\n');
    
    fbalt_upr = input('?');
    
    if (fbalt_upr >= fbalt_lwr)
        
        break;
        
    end
    
end

% request type of optimization

while(1)
    
    fprintf('\n\n      optimization menu\n');
    
    fprintf('\n <1> minimize departure delta-v\n');
    
    fprintf('\n <2> minimize arrival delta-v\n');
    
    fprintf('\n <3> minimize total delta-v\n');
    
    fprintf('\n selection (1, 2 or 3)\n');
    
    otype = input('?');
    
    if (otype == 1 || otype == 2 || otype == 3)
        
        break;
        
    end
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve the patched-conic gravity assist problem %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; home;

xg = zeros(3, 1);

% control variables initial guesses
% (departure, flyby, and arrival tdb julian dates)

xg(1) = jdtdbi1 - jdtdb0;

xg(2) = jdtdbi2 - jdtdb0;

xg(3) = jdtdbi3 - jdtdb0;

% bounds on control variables

xlwr(1) = xg(1) - ddays1;
xupr(1) = xg(1) + ddays1;

xlwr(2) = xg(2) - ddays2;
xupr(2) = xg(2) + ddays2;

xlwr(3) = xg(3) - ddays3;
xupr(3) = xg(3) + ddays3;

xlwr = xlwr';
xupr = xupr';

% bounds on objective function

flow(1) = -inf;
fupp(1) = +inf;

% bounds on flyby altitude

flow(2) = fbalt_lwr;
fupp(2) = fbalt_upr;

% bounds on v-infinity matching (equality) constraint

flow(3) = 0.0;
fupp(3) = 0.0;

flow = flow';
fupp = fupp';

xmul = zeros(3, 1);

xstate = zeros(3, 1);

fmul = zeros(3, 1);

fstate = zeros(3, 1);

% use snopt to find optimal solution

[x, f, inform, xmul, fmul] = snopt(xg, xlwr, xupr, xmul, xstate, ...
    flow, fupp, fmul, fstate, 'fbyfunc');

% evaluate optimal solution

[~, g] = fbyfunc(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute orientation of the departure hyperbola
% (Earth mean equator and equinox of J2000)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dveq1 = eq2000 * dv_launch;

dveqm1 = norm(dveq1);

decl1 = 90.0 - rtd * acos(dveq1(3) / dveqm1);

rasc1 = rtd * atan3(dveq1(2), dveq1(1));

% compute tdb julian dates of optimal transfer

jdtdb1 =  jdtdb0 + x(1);

jdtdb2 =  jdtdb0 + x(2);

jdtdb3 =  jdtdb0 + x(3);

% convert solution julian dates to calendar dates and tdb times

[cdstr1, tdbstr1] = jd2str(jdtdb1);

[cdstr2, tdbstr2] = jd2str(jdtdb2);

[cdstr3, tdbstr3] = jd2str(jdtdb3);

%%%%%%%%%%%%%%%%%
% print results %
%%%%%%%%%%%%%%%%%

fprintf('\nprogram flyby_matlab');

fprintf('\n--------------------\n');

fprintf('\npatched-conic gravity assist trajectory analysis\n');

if (otype == 1)
    
    fprintf('\nminimize departure delta-v\n');
    
end

if (otype == 2)
    
    fprintf('\nminimize arrival delta-v\n');
    
end

if (otype == 3)
    
    fprintf('\nminimize total delta-v\n');
    
end

fprintf('\ndeparture planet ');

disp(pname(ip1));

fprintf('flyby planet     ');

disp(pname(ip2));

fprintf('arrival planet   ');

disp(pname(ip3));

% departure conditions

fprintf('\nLAUNCH CONDITIONS');
fprintf('\n=================\n');

fprintf('\ndeparture calendar date        ');

disp(cdstr1);

fprintf('\ndeparture TDB time             ');

disp(tdbstr1);

fprintf('\ndeparture TDB julian date    %12.6f', jdtdb1);

fprintf('\n\nheliocentric launch delta-v vector and magnitude');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

fprintf('\nlaunch delta-vx            %12.6f meters/second', 1000.0 * dv_launch(1));
fprintf('\nlaunch delta-vy            %12.6f meters/second', 1000.0 * dv_launch(2));
fprintf('\nlaunch delta-vz            %12.6f meters/second', 1000.0 * dv_launch(3));

fprintf('\n\nlaunch delta-v             %12.6f meters/second\n', 1000.0 * dvm_launch);

fprintf('\nlaunch hyperbola characteristics');
fprintf('\n(Earth mean equator and equinox of J2000)');
fprintf('\n-----------------------------------------\n');

fprintf('\nasymptote right ascension  %12.6f degrees\n', rasc1);

fprintf('\nasymptote declination      %12.6f degrees\n', decl1);

fprintf('\nlaunch energy              %12.6f kilometer^2/second^2\n', c3launch);

fprintf('\npost-impulse heliocentric orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

oev1 = eci2orb1(smu, rito1, vito1);

oeprint1(smu, oev1, 3);

svprint(rito1, vito1);

fprintf('heliocentric orbital elements and state vector of the departure planet');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

[rp1, vp1] = p2000(ip1, jdtdb1);

oev1 = eci2orb1(smu, rp1, vp1);

oeprint1(smu, oev1, 3);

svprint(rp1, vp1);

% flyby conditions

fprintf('\nPATCHED-CONIC FLYBY CONDITIONS');
fprintf('\n==============================\n');

fprintf('\nflyby calendar date            ');

disp(cdstr2);

fprintf('\nflyby TDB time                 ');

disp(tdbstr2);

fprintf('\nflyby TDB julian date        %12.6f\n', jdtdb2);

fprintf('\nlaunch-to-flyby time       %12.6f days\n', jdtdb2 - jdtdb1);

fprintf('\nv-infinity in              %12.6f meters/second', 1000.0 * vinfm_in);

fprintf('\nv-infinity out             %12.6f meters/second\n', 1000.0 * vinfm_out);

fprintf('\nflyby altitude             %12.6f kilometers\n', fbalt);

% turn angles

tmp = rep(ip2) * vinfm_in * vinfm_in / pmu(ip2);

tangle1 = 2.0 * asin(1.0 / (1.0 + tmp));

fprintf('\nmaximum turn angle         %12.6f degrees', tangle1 * rtd);

tmp = (fbalt + rep(ip2)) * vinfm_in * vinfm_in / pmu(ip2);

tangle2 = 2.0 * asin(1.0 / (1.0 + tmp));

fprintf('\nactual turn angle          %12.6f degrees\n', tangle2 * rtd);

fprintf('\nheliocentric delta-v       %12.6f meters/second\n', 1000.0 * dvh);

% heliocentric deltavs

dvh_max = sqrt(pmu(ip2) / rep(ip2));

fprintf('\nmax heliocentric delta-v   %12.6f meters/second\n', 1000.0 * dvh_max);

% planet-centered orbital elements at periapsis

oev_periapsis = fbhyper(pmu(ip2), vinf_in, vinf_out, fbrp);

[rp2sc1, vp2sc1] = orb2eci(pmu(ip2), oev_periapsis);

fprintf('\nplanet-centered orbital elements and state vector of the spacecraft at periapsis');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

oeprint1(pmu(ip2), oev_periapsis, 1);

svprint(rp2sc1, vp2sc1);

[bplane, ~, ~, ~] = rv2bplane(pmu(ip2), rp2sc1, vp2sc1);

fprintf('b-plane coordinates of the spacecraft at periapsis');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

bpprint(rp2sc1, vp2sc1, bplane);

cdecl_asy = cos(bplane(2));
crasc_asy = cos(bplane(3));

srasc_asy = sin(bplane(3));
sdecl_asy = sin(bplane(2));

fprintf('\nheliocentric orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

tau = 86400.0 * (jdtdb2 - jdtdb1);

[rp, vp] = twobody2(smu, tau, rito1, vito1);

oev = eci2orb1(smu, rp, vp);

oeprint1(smu, oev, 3);

svprint(rp, vp);

fprintf('\nheliocentric orbital elements and state vector of the flyby planet');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

[rp, vp] = p2000(ip2, jdtdb2);

oev = eci2orb1(smu, rp, vp);

oeprint1(smu, oev, 3);

svprint(rp, vp);

% arrival conditions

fprintf('\nARRIVAL CONDITIONS');
fprintf('\n==================\n');

fprintf('\narrival calendar date          ');

disp(cdstr3);

fprintf('\narrival TDB time               ');

disp(tdbstr3);

fprintf('\narrival TDB julian date      %12.6f\n', jdtdb3);

fprintf('\nflyby-to-arrival time      %12.6f days\n', jdtdb3 - jdtdb2);

fprintf('\nheliocentric arrival delta-v vector and magnitude');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

fprintf('\narrival delta-vx            %12.6f meters/second', 1000.0 * dv_arrival(1));
fprintf('\narrival delta-vy            %12.6f meters/second', 1000.0 * dv_arrival(2));
fprintf('\narrival delta-vz            %12.6f meters/second', 1000.0 * dv_arrival(3));

fprintf('\n\narrival delta-v             %12.6f meters/second\n', 1000.0 * dvm_arrival);

% two-body propagation of the second leg of the transfer orbit

[rf3, vf3] = twobody2 (smu, 86400 * (jdtdb3 - jdtdb2), rito2, vito2);

fprintf('\npre-impulse heliocentric orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

oev3 = eci2orb1(smu, rf3, vf3);

oeprint1(smu, oev3, 3);

svprint(rf3, vf3);

fprintf('\npost-impulse heliocentric orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

rf4 = rf3;

vf4(1) = vf3(1) + dv_arrival(1);
vf4(2) = vf3(2) + dv_arrival(2);
vf4(3) = vf3(3) + dv_arrival(3);

oev4 = eci2orb1(smu, rf4, vf4);

oeprint1(smu, oev4, 3);

svprint(rf4, vf4);

fprintf('heliocentric orbital elements and state vector of the arrival planet');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

[rp, vp] = p2000(ip3, jdtdb3);

oev = eci2orb1(smu, rp, vp);

oeprint1(smu, oev, 3);

svprint(rp, vp);

fprintf('\nMISSION SUMMARY');
fprintf('\n===============\n');

fprintf('\ntotal delta-v              %12.6f meters/second\n', ...
    1000.0 * (dvm_launch + dvm_arrival));

fprintf('\ntotal energy               %12.6f kilometer^2/second^2\n', dvm_launch^2 + dvm_arrival^2);

fprintf('\ntotal mission duration     %12.6f days\n\n', jdtdb3 - jdtdb1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% two-body integration of the trajectory from
% first impulse to SOI entrance at flyby planet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set up for ode45

options = odeset('RelTol', 1.0e-6, 'AbsTol', 1.0e-8, 'Events', @soi_event);

% define maximum search time (seconds)

tof = 86400.0 * (jdtdb2 - jdtdb1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for flyby planet SOI entrance conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[~, ~, tevent, ~, ~] = ode45(@twobody_heqm, [0 tof], [rito1 vito1], options);

jdtdb_soi = jdtdb1 + tevent / 86400.0;

[cdstr_ca, utstr_ca] = jd2str(jdtdb_soi);

fprintf('\nPATCHED-CONIC SOI ENTRANCE CONDITIONS');
fprintf('\n=====================================\n');

fprintf('\ncalendar date          ');

disp(cdstr_ca);

fprintf('\nTDB time               ');

disp(utstr_ca);

fprintf('\nTDB julian date      %12.6f\n', jdtdb_soi);

fprintf('\nplanet-centered orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

oev_soi = eci2orb1(pmu(ip2), rp2sc, vp2sc);

oeprint1(pmu(ip2), oev_soi, 1);

svprint(rp2sc, vp2sc);

[bplane, ~, ~, ~] = rv2bplane(pmu(ip2), rp2sc, vp2sc);

fprintf('b-plane coordinates at sphere-of-influence');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

bpprint(rp2sc, vp2sc, bplane);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% integrate trajectory from SOI entrance
% to closest approach to the flyby planet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set up for ode45

options = odeset('RelTol', 1.0e-10, 'AbsTol', 1.0e-10, 'Events', @fpa_event);

% define maximum search time (seconds)

tof = 10.0 * 86400.0;

[~, ~, tevent, yevent, ie] = ode45(@peqm, [0 tof], [rp2sc vp2sc], options);

% extract state vector at closest approach

rsc_ca = yevent(1:3);

vsc_ca = yevent(4:6);

% B-plane coordinates of flyby hyperbola at closest approach

[bplane, rv, tv, ibperr] = rv2bplane(pmu(ip2), rsc_ca, vsc_ca);

% tdb julian date at closest approach

jdtdb_ca = jdtdb_soi + tevent / 86400.0;

[cdstr_ca, utstr_ca] = jd2str(jdtdb_ca);

fprintf('\n\nNUMERICALLY INTEGRATED CLOSEST APPROACH CONDITIONS');
fprintf('\n==================================================\n');

fprintf('\ncalendar date          ');

disp(cdstr_ca);

fprintf('\nTDB time               ');

disp(utstr_ca);

fprintf('\nTDB julian date      %12.6f\n', jdtdb_ca);

fprintf('\nplanet-centered orbital elements and state vector of the spacecraft');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

oev_ca = eci2orb1(pmu(ip2), rsc_ca, vsc_ca);

oeprint1(pmu(ip2), oev_ca, 1);

svprint(rsc_ca, vsc_ca);

fprintf('b-plane coordinates at closest approach');
fprintf('\n(Earth mean ecliptic and equinox of J2000)');
fprintf('\n------------------------------------------\n');

bpprint(rsc_ca, vsc_ca, bplane);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% heliocentric and flyby graphics %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

while(1)
    
    fprintf('\n\nwould you like to create graphics for this mission (y = yes, n = no)\n');
    
    slct = input('?', 's');
    
    if (slct == 'y' || slct == 'n')
        break;
    end
    
end

if (slct == 'y')
    
    while(1)
        
        fprintf('\n\nplease input the plot step size (days)\n');
        
        deltat = input('?');
        
        if (deltat > 0)
            break;
        end
        
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plot heliocentric planetary orbits and transfer trajectory
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    figure(1);
    
    % departure planet semimajor axis and period
    
    [r1, v1] = p2000(ip1, jdtdb1);
    
    oev1 = eci2orb1(smu, r1, v1);
    
    period1 = 2.0 * pi * oev1(1) * sqrt(oev1(1) / smu) / 86400.0;
    
    % flyby planet semimajor axis and period
    
    [r2, v2] = p2000(ip2, jdtdb1);
    
    oev2 = eci2orb1(smu, r2, v2);
    
    period2 = 2.0 * pi * oev2(1) * sqrt(oev2(1) / smu) / 86400.0;
    
    % arrival planet semimajor axis and period
    
    [r3, v3] = p2000(ip3, jdtdb1);
    
    oev3 = eci2orb1(smu, r3, v3);
    
    period3 = 2.0 * pi * oev3(1) * sqrt(oev3(1) / smu) / 86400.0;
    
    % number of "planet" points to plot
    
    npts1 = fix(period1 / deltat);
    
    npts2 = fix(period2 / deltat);
    
    npts3 = fix(period3 / deltat);
    
    % departure planet orbit
    
    for i = 0:1:npts1
        
        jdtdb = jdtdb1 + i * deltat;
        
        [r1, ~] = p2000(ip1, jdtdb);
        
        x1(i + 1) = r1(1) / aunit;
        
        y1(i + 1) = r1(2) / aunit;
        
    end
    
    % compute last data point
    
    [r1, v1] = p2000(ip1, jdtdb1 + period1);
    
    x1(npts1 + 1) = r1(1) / aunit;
    
    y1(npts1 + 1) = r1(2) / aunit;
    
    % flyby planet heliocentric orbit
    
    for i = 0:1:npts2
        
        jdtdb = jdtdb1 + i * deltat;
        
        [r2, ~] = p2000(ip2, jdtdb);
        
        x2(i + 1) = r2(1) / aunit;
        
        y2(i + 1) = r2(2) / aunit;
        
    end
    
    % compute last data point
    
    [r2, v2] = p2000(ip2, jdtdb1 + period2);
    
    x2(npts2 + 1) = r2(1) / aunit;
    
    y2(npts2 + 1) = r2(2) / aunit;
    
    % arrival planet heliocentric orbit
    
    for i = 0:1:npts3
        
        jdtdb = jdtdb1 + i * deltat;
        
        [r3, ~] = p2000(ip3, jdtdb);
        
        x3(i + 1) = r3(1) / aunit;
        
        y3(i + 1) = r3(2) / aunit;
        
    end
    
    % compute last data point
    
    [r3, v3] = p2000(ip3, jdtdb1 + period3);
    
    x3(npts3 + 1) = r3(1) / aunit;
    
    y3(npts3 + 1) = r3(2) / aunit;
    
    % first leg of transfer orbit
    
    npts4 = fix((jdtdb2 - jdtdb1) / deltat);
    
    for i = 0:1:npts4
        
        tau = 86400.0 * i * deltat;
        
        [rft1, ~] = twobody2(smu, tau, rito1, vito1);
        
        x4(i + 1) = rft1(1) / aunit;
        
        y4(i + 1) = rft1(2) / aunit;
        
    end
    
    % compute last data point of first leg
    
    tau = 86400.0 * (jdtdb2 - jdtdb1);
    
    [rft1, vft1] = twobody2(smu, tau, rito1, vito1);
    
    x4(npts4 + 1) = rft1(1) / aunit;
    
    y4(npts4 + 1) = rft1(2) / aunit;
    
    % second leg of heliocentric transfer orbit
    
    npts5 = fix((jdtdb3 - jdtdb2) / deltat);
    
    for i = 0:1:npts5
        
        tau = 86400.0 * i * deltat;
        
        [rft2, vft2] = twobody2(smu, tau, rito2, vito2);
        
        x5(i + 1) = rft2(1) / aunit;
        
        y5(i + 1) = rft2(2) / aunit;
        
    end
    
    % compute last data point of second leg
    
    tau = 86400.0 * (jdtdb3 - jdtdb2);
    
    [rfto2, vfto2] = twobody2(smu, tau, rito2, vito2);
    
    x5(npts5 + 1) = rfto2(1) / aunit;
    
    y5(npts5 + 1) = rfto2(2) / aunit;
    
    % determine vernal equinox "size"
    
    xve = oev1(1) / aunit;
    
    if (oev2(1) > oev1(1))
        
        xve = oev2(1) / aunit;
        
    end
    
    if (oev3(1) > oev2(1))
        
        xve = oev3(1) / aunit;
        
    end
    
    % plot heliocentric planet orbits and patched-conic transfer trajectory
    
    hold on;
    
    plot(x1, y1, '.b');
    
    plot(x1, y1, '-b');
    
    plot(x2, y2, '.g');
    
    plot(x2, y2, '-g');
    
    plot(x3, y3, '.r');
    
    plot(x3, y3, '-r');
    
    plot(x4, y4, '.k');
    
    plot(x4, y4, '-k');
    
    plot(x5, y5, '.k');
    
    plot(x5, y5, '-k');
    
    plot(0, 0, 'hy', 'MarkerSize', 10);
    
    % plot and label vernal equinox direction
    
    line ([0.05, xve], [0, 0], 'Color', 'black');
    
    text(1.05 * xve, 0, '\Upsilon');
    
    % label launch, flyby and arrival locations
    
    plot(x4(1), y4(1), '*k');
    
    text(x4(1) + 0.06, y4(1) - 0.01, 'D');
    
    plot(x5(1), y5(1), '*k');
    
    text(x5(1) + 0.06, y5(1) + 0.01, 'F');
    
    plot(x5(npts5 + 1), y5(npts5 + 1), '*k');
    
    text(x5(npts5 + 1) + 0.06, y5(npts5 + 1) - 0.01, 'A');
    
    axis equal;
    
    % display launch, flyby and arrival dates
    
    text(0.85 * xve, -xve + 0.8, 'Departure', 'FontSize', 8);
    
    text(0.875 * xve, -xve + 0.7, pname(ip1), 'FontSize', 8);
    
    text(0.875 * xve, -xve + 0.6, cdstr1, 'FontSize', 8);
    
    text(0.85 * xve, -xve + 0.4, 'Arrival', 'FontSize', 8);
    
    text(0.875 * xve, -xve + 0.3, pname(ip3), 'FontSize', 8);
    
    text(0.875 * xve, -xve + 0.2, cdstr3, 'FontSize', 8);
    
    text(0.85 * xve, xve - 0.2, 'Flyby', 'FontSize', 8);
    
    text(0.875 * xve, xve - 0.3, pname(ip2), 'FontSize', 8);
    
    text(0.875 * xve, xve - 0.4, cdstr2, 'FontSize', 8);
    
    text(0.875 * xve, xve - 0.5, tdbstr2, 'FontSize', 8);
    
    % label plot and axes in astronomical units (AU)
    
    xlabel('X coordinate (AU)', 'FontSize', 12);
    
    ylabel('Y coordinate (AU)', 'FontSize', 12);
    
    title('Patched-conic Gravity Assist Trajectory', 'FontSize', 16);
    
    zoom on;
    
    % the next line creates a color eps graphics file with tiff preview
    
    print -depsc -tiff -r300 flyby_matlab1.eps
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % create three-dimensional graphics of the flyby
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % unit pointing vector from the flyby planet to the sun
    
    [rp, vp] = p2000(ip2, jdtdb2);
    
    up2s(1) = -rp(1) / norm(rp);
    up2s(2) = -rp(2) / norm(rp);
    up2s(3) = -rp(3) / norm(rp);
    
    sun_axisx = [0 5 * up2s(1)];
    sun_axisy = [0 5 * up2s(2)];
    sun_axisz = [0 5 * up2s(3)];
    
    % unit velocity vector of the flyby planet
    
    uv(1) = vp(1) / norm(vp);
    uv(2) = vp(2) / norm(vp);
    uv(3) = vp(3) / norm(vp);
    
    uv_axisx = [0 5 * uv(1)];
    uv_axisy = [0 5 * uv(2)];
    uv_axisz = [0 5 * uv(3)];
    
    dtof = 5000.0;
    
    [ri_ho, vi_ho] = twobody2 (pmu(ip2), -dtof, rp2sc1, vp2sc1);
    
    deltat1 = 2.0 * dtof / 300;
    
    simtime1 = -deltat1;
    
    for i = 1:1:301
        
        simtime1 = simtime1 + deltat1;
        
        % hyperbolic orbit "normalized" position vector
        
        [rf, vf] = twobody2 (pmu(ip2), simtime1, ri_ho, vi_ho);
        
        rp1_x(i) = rf(1) / rep(ip2);
        
        rp1_y(i) = rf(2) / rep(ip2);
        
        rp1_z(i) = rf(3) / rep(ip2);
        
    end
    
    % create axes vectors
    
    xaxisx = [1 1.5];
    xaxisy = [0 0];
    xaxisz = [0 0];
    
    yaxisx = [0 0];
    yaxisy = [1 1.5];
    yaxisz = [0 0];
    
    zaxisx = [0 0];
    zaxisy = [0 0];
    zaxisz = [1 1.5];
    
    figure(2);
    
    hold on;
    
    grid on;
    
    % plot flyby planet
    
    [x, y, z] = sphere(24);
    
    h = surf(x, y, z);
    
    colormap([127/255 1 222/255]);
    
    set (h, 'edgecolor', [1 1 1]);
    
    % plot coordinate system axes
    
    plot3(xaxisx, xaxisy, xaxisz, '-r', 'LineWidth', 1);
    
    plot3(yaxisx, yaxisy, yaxisz, '-g', 'LineWidth', 1);
    
    plot3(zaxisx, zaxisy, zaxisz, '-b', 'LineWidth', 1);
    
    % plot unit pointing vector to the sun
    
    plot3(sun_axisx, sun_axisy, sun_axisz, '-y', 'LineWidth', 1);
    
    plot3(uv_axisx, uv_axisy, uv_axisz, '-k', 'LineWidth', 1);
    
    % plot planet-centered flyby hyperbolic orbit
    
    plot3(rp1_x, rp1_y, rp1_z, '-m', 'LineWidth', 1);
    
    % label incoming leg of flyby hyperbola
    
    plot3(rp1_x(1), rp1_y(1), rp1_z(1), '*m');
    
    % label periapsis of flyby hyperbola
    
    plot3(rp2sc1(1) / rep(ip2), rp2sc1(2) / rep(ip2), rp2sc1(3) / rep(ip2), 'om');
    
    % label plot in flyby planet radii (PR)
    
    xlabel('X coordinate (PR)', 'FontSize', 12);
    
    ylabel('Y coordinate (PR)', 'FontSize', 12);
    
    zlabel('Z coordinate (PR)', 'FontSize', 12);
    
    title('Patched-conic Gravity Assist Trajectory', 'FontSize', 16);
    
    axis equal;
    
    view(38, 30);
    
    rotate3d on;
    
    % the next line creates a color eps graphics file with tiff preview
    
    print -depsc -tiff -r300 flyby_matlab2.eps
    
end

Error:

Index exceeds the number of array elements (0).

Error in Untitled2 (line 742)

rsc_ca = yevent(1:3);

Answer 1

[Image Analyst]

yevent is empty - it has no elements. Debug your program to find out why.

Debugging in MATLAB | Doug's MATLAB Video Tutorials

[Walter Roberson]

In particular, yevent will be empty if the call runs to the last time given, or if the call stops early because it cannot integrate (function is too steep or has a discontinuity.)

You did not happen to include peqm so we cannot test.