Found a code that works!
clc
clear all
% Create network for curve fitting
hiddenLayerSize = 4; % Number of intermediate network neurons
net = fitnet(hiddenLayerSize);
WB = getwb(net); % Only net.b{1} = zeros(10,1)is defined
rng(4151945); % Initialize the RNG so that results can be duplicated
M = [1:1:10];
M = [M,M,M,M,M].*rand();
M = [M,M].*rand();
M = [M,M,M,M,M].*10;
M = [M,M].*10;
a=M.*rand().*2^rand()+5*rand()-5*rand();
b=M.*rand().*2^rand()+5*rand()-5*rand();
c=M.*rand().*2^rand()+5*rand()-5*rand();
n=rand(1,1000)*0.05;
y = 5*a + b.*c + 7*c + n;
x=[a; b; c];
t=y;
% Setting the sample size
net.divideFcn = 'dividerand'; % Split random data
net.divideMode = 'sample';
net.divideParam.trainRatio = 70/100;
net.divideParam.valRatio = 15/100;
net.divideParam.testRatio = 15/100;
net = train(net,x,t);
[xn, xsettings] = mapminmax(x);
[tn, tsettings] = mapminmax(t);
% syms p q r real
% X = [p,q,r]';
X = [22,25,21]'
x000 = mapminmax('apply',X,xsettings);
b1 = net.b{1};
b2 = net.b{2};
IW = net.IW{1,1};
LW = net.LW{2,1};
yn = b2 + LW * tanh(b1 + (IW * x000))
y2 = mapminmax.reverse(yn,tsettings)
y1 = sim(net,X)
y1compare = 5*X(1) + X(2)*X(3) + 7*X(3)