Gauss-Legendre Matlab: How to use Newton method to approximate k_1 and k

Recently as a project I have been working on a program in matlab designed to implement the Gauss-Legendre order 4 method of solving an IVP. Numerical analysis- and coding in particular- is somewhat of a weakness of mine, and thus it has been rather tough going. I used the explicit Euler method to initially seed the method, but I've been stuck for a very long time on how to use the Newton method to get closer values of k1 and k2.

Thusfar, my code looks as follows:

%Gauss-Butcher Order 4
function[y] = GBOF(f,fprime,y_0,maxit,ertol)
A = [1/4,1/4-sqrt(3)/6;1/4 sqrt(3)/6,1/4];
h = [1,0,0,0,0,0,0,0,0,0,0,0,0];
t = [0,0,0,0,0,0,0,0,0,0,0,0,0];
for n = 1:12
    h(n 1) = 1/(2^n);
    t(n 1) = t(n) h(n);
end
y = zeros(size(t));
y(1) = y_0;
niter = 1;