Is there an input/output relationship among the 2 variables? That is, does changing one affects the other in a cause-effect manner? If so, you need to separate out the variables into input data and output data. If not, you can treat the data set as 3D time series. In either case, start by packaging the data as an IDDATA object. In the following I assume you have the time series case:
z = iddata(data, [], 1); % data has 3 columns and 50 rows
% compute linear ARX model of some order NA
NA = [3 0 0; 0 4 0; 0 0 5];
m1 = arx(z, NA);
% compute a nonlinear model that uses tree partition nonlinearity and same orders
m2 = nlarx(z, NA, 'treepartition');
% compute a nonlinear model that employs sigmoid network with 5 units for all three time series:
m3 = nlarx(z, NA, sigmoidnet('num', 3))
Using diagonal NA means the three variables are treated as independent (not influencing each other). Introduce nonzero values for the nondiagonal entries of NA to realize a coupled model. NA values would be determined using prior information (knowledge of what the timeseries represent and how they might be influencing each other) or numerical order search (try various combinations and see which one(s) give good fit to an independent validation dataset).