Fit two datasets with shared and individual parameters to a function containing a numerical integration

In this particular problem I want to write some code that fits function 'fun1' to two different datasets (data1 and data2). My fitting function 'fun1' comes from a numerical integration and the independent variable 'X' is changes the integration boundaries. Another help post I cannot recall suggested the implementation with arrayfun which worked fine for me. When all fitting parameters (p1 to p4) were independent I could find a fitting function for 'data1' and 'data2' separately. I suspect that in my probem only p1 differs between the two datasets, so I was trying to rewrite the code according to this post. I recognise that I need to add another variable to arrayfun, since p1 is a vector now, but I am getting an error that 'all of the input arguments must be of the same size and shape'.

X=[data1{:,1}/max(data1{:,1}); data2{:,1}/max(data2{:,1})]
Y=[data1{:,2}/max(data1{:,2}); data2{:,2}/max(data2{:,2})]
dsid=[1*ones(length(X(:,1))/2,1); 2*ones(length(X(:,1))/2,1)]
gscatter(X,Y,dsid)

fun= @(p1,p2,p3,p4,r) arrayfun(@(dvec,R) p4+p2*integral(@(x) exp(-p3*x).*erf(p3.*(sqrt(R.^2-(p1(dvec)-x).^2)))/(2*p3),p1(dvec)-R,p1(dvec)+R),dsid,r);
fun1= @(params,X) fun(params(1:2),params(3),params(4),params(5),X);

x01=[5,12,0.3,0.3,0.1];
options = optimoptions(@lsqcurvefit,'MaxFunctionEvaluations',100000,'MaxIterations',20000,'Algorithm','trust-region-reflective');
params = lsqcurvefit(fun1, x01, X, Y,[],[],options);