Gaussian fit fails converging to a good result

[LStarM]

Hi, my name is Lorenzo. I am not sure this is the right place to ask for support but I try. I am working on a project for my research where I have to perform a large number of fit on different images. Matlab fitting routines are too slow for this task also running in parallel. So searching online for help I found gpufit. I would like to fit my own model but as a test I started with a simple gaussian fit but I am struggling to get it working. I do not get any error but the fit does not converge at a reasonable result even if the starting point is very close to the expected values. Also those curves for which more iterations are done remain very far from the expected result, sometimes it does not even move from the initial guess. I have also tried starting further from a good initial guess to avoid being in a local minimum but still not getting any good result. Not sure what the problem is but it does not seem to be math related.
I attach the full working Matlab code. Can anyone help?
I am working with CUDA 11.2 and Matlab 2022b.

g = gpuDevice(1);
reset(g);

% generate x axis
x = -1:0.02:1;
xfine = x(1):0.0001:x(end);

% generate gaussians
data = 200exp(-(x-0.2).^2/(20.1^2))+10;
data = repmat(data,[1e4,1]);
data = normrnd(data,sqrt(data));

% generate initial guess
p0 = max(data,[],2)-min(data,[],2);
p1 = sum(data.x,2)./sum(data,2);
p2 = (abs(sum((x-p1).^2.(data),2)./sum(data,2)));
p3 = min(data,[],2);

% plot starting point
idx = 500;
subplot(1,2,1);
plot(x,data(idx,:)); hold on; plot(xfine,p0(idx)exp(-(xfine-p1(idx)).^2/(2p2(idx)^2))+p3(idx)); hold off

% arrange data and initial guess in the required format
data = single(data)';
weigths = [];
initialParam = single([p0';p1';p2';p3']);
model_id = 0;
estimator_id = 0;
maxNiterations = 25;
tolerance = 1e-4;

% perform fit
[parameters, states, chi_squares, n_iterations, time] = gpufit(data, [], model_id, initialParam, tolerance, maxNiterations, [], estimator_id, x);

% plot
subplot(1,2,2)
plot(x,data(:,idx)); hold on; plot(xfine,parameters(1,idx)exp(-(xfine-parameters(2,idx)).^2/(2parameters(3,idx)^2))+parameters(4,idx)); hold off

reset(g);
clear GpufitMex;