Something wrong with my implementation of the bias/variance diagnostic in polynomial regression I'm trying to diagnosing bias/variance so I have the below Octavecode:
%========== Begin - constants declaration ==========% 
x_training_percent = 0.7;
cv_set_percent = 0.3;
%========== End - constants declaration ==========% 

load('data/dataset.m');
data = [X y];
data = data(:, randperm(size(data,2)));


[m, n] = size(X);% m: Number of examples, n: Number of features.

%========== Begin - Getting traingin and CV sets ==========% 
training_set_size = round(m * x_training_percent);
cv_set_size = round(m * cv_set_percent);
%test_set_size = round(m * test_set_percent);

x_training_o = data(1:training_set_size, 1:(end - 1));
y_training = data(1:training_set_size, end);

x_cv_o = data((training_set_size + 1):(training_set_size + cv_set_size), 1:(end - 1));
y_cv = data((training_set_size + 1):(training_set_size + cv_set_size), end);

%========== End - Getting traingin and CV sets ==========% 

max_p = 20; % Max degree polynomial

cv_error = zeros(max_p, 1);
training_error = zeros(max_p, 1);

for p = 1:max_p

  % Processing training set
  x_training = x_training_o;
  x_training = polyFeatures(x_training, p); % Adding polynomial terms from 1 to p
  x_training = [ones(training_set_size, 1) x_training];

  % Processing cross validation set
  x_cv = x_cv_o;
  x_cv = polyFeatures(x_cv, p); % Adding polynomial terms from 1 to p
  x_cv = [ones(cv_set_size, 1) x_cv];

  %========== Begin - Training ==========% 
  lambda = 0
  theta = trainLinearReg(x_training, y_training, lambda);
  %========== End - Training ==========% 

  %========== Begin - Computing prediction errors with polinomial degree p ==========% 
    predictions = x_training * theta; % Predictions with training set
    training_error(p, :) = (1 / (2 * training_set_size)) * sum((predictions - y_training) .^ 2);

    cv_predictions = x_cv * theta; % Predictions with cross validation set
    cv_error(p, :) = (1 / (2 * cv_set_size)) * sum((cv_predictions - y_cv) .^ 2);
  %========== End - Computing prediction errors ==========% 


end
plot(1:max_p, training_error, 1:max_p, cv_error);

legend('Train', 'Cross Validation')
xlabel('Degree of Polynomial')
ylabel('Error')

My outputs is:

But it's very different from the following result shown by Andrew Ng:

So what am I doing wrong here?
P.D: You can see the complete source code here
 A: Here is how you should debug your code:


*

*For now, forget about cross-validation. Your most serious problem is the training error increasing with the polynomial degree, while it should decrease (and become zero for big enough polynomial degree). Make the code without any traces of cross-validation. It will be much easier to check.

*If the error will not be reproduced (when you remove the cross-validation from your code, the training error suddenly decreases), then add the cross-validation part gradually (small chunks of code) until you find which part spoils the correct work, and understand why. Otherwise, go to step 3.

*Now you have a code without cross-validation, which does not work properly. Make the simulated data which ideally fits, say, quadratic polynomial, and check if it works properly. If it doesn't, try the linear polynomial. At the end you will have something very simple which does not work.  You will probably be able to figure out why.

*Posting your whole code which does not work and asking people "what did I do wrong" is a bad idea. It means that you want people to debug your code instead of you.  Why should they do it?  You should first try to make your code as small as possible while still retaining the error.  If you can remove something from your code, and the error is still here, then the code is NOT as simple as possible.  If there is a part of your code which you did not try to remove, the code is NOT as simple as possible.
A: I had a very similar problem before, and it was because the step size in my gradient descent implementation was too large and the tolerance too forgiving.
As you add higher order features their magnitude may be rapidly increasing, making deviations from the true coefficient values cause greater error. The fact that your training error increases as you add higher order features suggests this may be happening.
A: I have added feature normalization (mean = 0 and standard deviation = 1) and seems to be that I found the problem. See the new curve (Just training set). However, error is so high 

And the code below...
addpath ./common;
addpath ./common/minFunc_2012/minFunc;
addpath ./common/minFunc_2012/minFunc/compiled;

load('data/dataset.m');
data = [X y];
data = data(:, randperm(size(data,2)));


[m, n] = size(X);% m: Number of examples, n: Number of features.

%========== Begin - Getting traingin and CV sets ==========% 
x_training_o = data(1:m, 1:(end - 1));
y_training = data(1:m, end);
%========== End - Getting traingin and CV sets ==========% 

max_p = 10; % Max degree polynomial

training_error = zeros(max_p, 1);

csvwrite("log/Xy_tr.txt", [x_training_o y_training]);

for p = 1:max_p
  x_training = x_training_o;
  x_training = polyFeatures(x_training, p); % Adding polynomial terms from 1 to p
  [x_training, mu, sigma] = featureNormalize(x_training);
  x_training = [ones(m, 1) x_training];


  %========== Begin - Training ==========% 
  lambda = 0
  theta = trainLinearReg(x_training, y_training, lambda);
  %========== End - Training ==========% 

  %========== Begin - Computing prediction errors with polinomial degree p ==========% 
    predictions = x_training * theta; % Predictions with training set

  training_error(p, :) = (1 / (2 * m)) * sum((predictions - y_training) .^ 2);

  %========== End - Computing prediction errors ==========% 

end

A: I suspect you may find better results when you use the ordinary least squares estimator instead of gradient descent:
$$\hat{ \beta } = (X'X)^{-1} X'Y.$$
