It is not clear to me exactly what type of model you hope to fit. Some people say 'multivariate regression' to mean there are multiple dependent variables. For examples, you might want to predict both humidity and temp from day. Other people use it to just mean one outcome but multiple predictors. For example, you could predict temperature from both day and humidity. Here is an instructional solution showing how you could get the parameter estimates, standard errors, and new predicted values from a multiple regression model predicting the outcome, temperature, from day and humidity (as well as the constant term).
# your original data
day = [4, 5, 6, 8]
temp = [97, 100, 98, 80]
humidity = [62, 46, 50, 55]
# create the design matrix
# intercept (1s), day and humidity as predictors
X = [1, 1, 1, 1; day; humidity]'
# linear parameter estimates
b = inv(X'*X)*X'*temp'
R = temp' - (X * b)
# residual variance
v = (R'*R)/(4 - 3)
# variance covariance matrix of parameters
Sigma = v * inv(X'*X)
# standard errors of parameters (b vector)
se = sqrt(diag(Sigma))
# new data for prediction, constant
# day is 7 and 9, humidity is 80
newdata = [1, 7, 80; 1, 9, 80]
# predicted values for day 7 and 9
pred = newdata * b
In practice, that is recreating the wheel, but since you are new to octave (and maybe regression?) I thought it might be helpful. Here is the simple way using built in functions to directly get the coefficients.
which would be
b from above, and you could postmultiply by new data (in your case day 7 and 9) to get predicted ("forecasts") values.