function y = fconc31(x)

global b1 b2 b3 b4
global c31
global cnew

b1 = x(1);
b2 = x(2);
b3 = x(3);
b4 = x(4);
[t,cnew] = ode45('freact1',[0:0.02:2],[1 0 0]);
c31new = cnew(round(linspace(1,101,31)),:);

y = sum(sum((c31new-c31).*(c31new-c31)));