How to aproximate a power demand code to reality?

Question

I've written a code for the usage of some cooking appliances, that will run at breakfast, lunch and dinner time for a 5 million population (1 million households). The aim is to plot the power demand in 24 hours.

This is the code i've made:

function [monthlyEnergyMWhcook,powerConsumptionMWcook] = Cooking()
% Constants for appliances (Watts)
ovenPower = 3000; % Watts (household appliance)
grillPower = 6000; % Watts (household appliance)
kettlePower = 2200; % Watts (household appliance)
toasterPower = 850; % Watts (household appliance)
microwavePower = 800; % Watts (household appliance)

% Daily usage durations in minutes, converted to number of intervals
ovenDuration = ceil((132 / 60) * (241 / 24)); % 2 hours and 12 minutes
grillDuration = ceil((126 / 60) * (241 / 24)); % 2 hours and 6 minutes
kettleDuration = ceil((11 / 60) * (241 / 24)); % 11 minutes
toasterDuration = ceil((9 / 60) * (241 / 24)); % 9 minutes
microwaveDuration = ceil((11 / 60) * (241 / 24)); % 11 minutes

% Time vector for 24-hour period (241 time steps)
timeVec = linspace(0, 24, 241);

% Initialize power consumption array
powerConsumption = zeros(1, length(timeVec));

% Define key meal preparation start times
breakfastStart = find(timeVec >= 7, 1); % 7 AM
lunchStart = find(timeVec >= 12, 1); % 12 PM
dinnerStart = find(timeVec >= 18, 1); % 6 PM

% Number of households
households = 1000000;  % Example number of households

% Appliance usage during typical meal times
% Breakfast usage
powerConsumption(breakfastStart:breakfastStart+kettleDuration-1) = (kettlePower+toasterPower+microwavePower) * households;
% Lunch usage
powerConsumption(lunchStart:lunchStart+ovenDuration-1) = ovenPower * households;

% Dinner usage
powerConsumption(dinnerStart:dinnerStart+grillDuration-1) = grillPower * households;

% Convert power consumption to MW
powerConsumptionMWcook = powerConsumption / 1000000; % Convert Watts to kW and then to MW

% Calculate total daily energy in MWh
totalEnergyMWh = sum(powerConsumptionMWcook) * (6 / 60); % Convert MW to MWh, accounting for 6 minutes per interval

% Print the total daily energy consumption
fprintf('Total daily energy consumed by cooking: %.2f MWh\n', totalEnergyMWh);
%Plotting the 24-hour power consumption for cooking
figure;
plot(timeVec, powerConsumptionMWcook);
xlabel('Time (Hours)');
ylabel('Power Consumption (MW)');
title('24-Hour Cooking Power Consumption in MW');
grid on;
%Constants
daysInMonth = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];  % Days in each month

% Daily energy consumption in MWh (from previous calculation)
dailyEnergyMWh = sum(powerConsumptionMWcook) * (6 / 60); % Recalculated to ensure accuracy

% Initialize monthly energy consumption array
monthlyEnergyMWhcook = zeros(1, 12);

% Calculate monthly energy consumption
for i = 1:12
    monthlyEnergyMWhcook(i) = dailyEnergyMWh * daysInMonth(i);
end

Problem:

The thing is that my code starts every appliance at the same time, it means for example the million kettles that will run for 12 mins start at the same time, so it generates huge peaks in power demand that do not represent reality.

Does anyone know how to make that during breakfast time for example between 8 and 9, all the kettles complete their 12 mins of usage, but starting at different time steps?

This is the graph obtained from the code

I am thinking that the graph should look like a gausssian curve but not sure.

I tried to model gaussian curves but the total diary energy increased quite a lot.

I've been searching how to perform my idea, but i haven't found a solution yet.

I would appreciate any kind of help, Kind regards.