Find proportionality using nlinfit

Hello good ppl!

I have an array x and an array y. y is presented below and is values of lumiosity as a function of days (x) for a supernova. I would like to find the proportionality of luminosity to time, L \prop x^(?). Here are my failed attempts using nlinfit. Also, y values are part of another array Lbolos.

 if true
x = 960:1236
y = Lbolos(261:537)
model = @(beta, s)( ( Lbolos(261)/960^(beta) )*s^(beta) ) % Lbolos(261)/960^(beta) is for normalisation/start at the right place
initial = [-4]
coeff = nlinfit(x,y, model, initial)
 end

if true
x = 960:1236
y = Lbolos(261:537)
model = @(beta, s)( beta(1)*s^(beta(2)) )
initial = [Lbolos(261)/(960^(-4)) -4]
coeff = nlinfit(x,y, model, initial)
end

both times the error message is: Error using nlinfit Error evaluating model function

Help is much appreciated, also a full example if possible<3

if true
  y = 1×277    
  7.436e+41   7.3707e+41   7.3046e+41   7.2376e+41     7.17e+41   7.1017e+41   7.0328e+41   6.9633e+41   6.8933e+41   6.8229e+41    6.752e+41   6.6808e+41   6.6093e+41   6.5374e+41   6.4654e+41   6.3931e+41   6.3207e+41   6.2482e+41   6.1757e+41   6.1031e+41   6.0305e+41    5.958e+41   5.8855e+41   5.8132e+41   5.7411e+41   5.6691e+41   5.5973e+41   5.5258e+41   5.4546e+41   5.3836e+41    5.313e+41   5.2427e+41   5.1729e+41   5.1034e+41   5.0344e+41   4.9658e+41   4.8977e+41   4.8301e+41   4.7631e+41   4.6965e+41   4.6305e+41   4.5651e+41   4.5002e+41    4.436e+41   4.3723e+41   4.3093e+41   4.2469e+41   4.1852e+41   4.1241e+41   4.0636e+41 ...

end