Hello I have to write a code for the following second order ODE: d^2/dx^2=7x-8x^2+2y. This is what I have so far but I keep getting- Attempted to access y(2); index out of bounds because numel(y)=1. I would appreciate some help

function dy=dydx(x,y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Create a system of fisrt order equations from the given second order which
%is given by the m.file dydx.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%let y1=y      and      y2=y'
%therefore: y1'=y'=y2        and       y2'=y''=2y1+8x(9-x)
dy=[y(2);2*y(1)+72*x-8*(x.^2)];

function [x,y]=eulermethod(dydx,xvec,y0,h)
%dydx is the name of the m file that evaluates the given second order ODE
%xvec=[xi;xf] where xi and xf are the initial and final values of the
%independent variable
%y0= initial value of the dependent variable
%h=step size
%[x,y]=output vector of independent and dependent variables respectively
xi=xvec(1);
xf=xvec(2);
x=(xi:h:xf);
n=length(x);
if x(n)<xf %If required an additional value of x is added so that the range is from xi to xf
    x(n+1)=xf;
    n=n+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Implementing Euler method 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y(1,:) = y0(:); % Start y at the initial value .
for i=2:n
    hh=x(i+1)-x(i);
    k=feval(dydx,x(i),y(:,i))';
    y(i+1,:)=y(i, :)+k*hh;
end

Answer 1

[James Tursa]

Here you got the states set up in rows:

y(1,:) = y0(:);

Yet here you have the call made as if the states are in columns:

    k=feval(dydx,x(i),y(:,i))';

So pick one or the other and be consistent.