When you use the Vectorized option for odeset, your ode function will be called with an unpredetermined number of columns for y, each corresponding to what would be one y vector if the function had been called without Vectorization turned on. The number of columns could vary between calls to your ode function.
The number of rows for y for your ode function will match length(y0). In your case your y0 is of length 2. But your code in f attempts to access right up to row ((2*N-1)+1) = 2*N . As your y0 is of length 2, your N would have to be 1 for that to work.
The number of rows you output must be the same as the number of rows for the input y. If you have 2*N rows of output you had better make your y0 length 2*N:
y0 = zeros(2*N,1);
Also, in order for your code to be able to work properly, your boundary conditions need to be consistent. You have f(0,t)=1, g(z,0)=0 and f(z,0)=0 so f(0,0)=1 due to the first condition but f(0,0)=0 due to the third condition. Your boundaries are therefore not consistent so after you fix your y0, you will not be able to achieve a useful output.
